Kode di Hide dalam default, untuk menampilkan
kode, klik Code .
Analisis sebelum ini : https://rpubs.com/ZenR_Prog/MPDW-Prak10
# -=( Install & Load Package Function )=-
install_load <- function (package1, ...) {
# convert arguments to vector
packages <- c(package1, ...)
# start loop to determine if each package is installed
for(package in packages){
# if package is installed locally, load
if(package %in% rownames(installed.packages()))
do.call('library', list(package))
# if package is not installed locally, download, then load
else {
install.packages(package)
do.call("library", list(package))
}
}
}
#Path Function
path <- function(){
gsub ( "\\\\", "/", readClipboard () )
}
#Copy path, Panggil function di console
#Copy r path, paste ke var yang diinginkan
#Export chart
export.chart <- "C:/Users/Fathan/Documents/Obsidian Vault/2. Kuliah/Smt 5/6. Metode Peramalan Deret Waktu/@Proj/STA1341-MPDW/Pertemuan 10-14/Chart2"
1 Pendahuluan
1.1 NVIDIA
Sumber : finance.yahoo.com
Terdapat beberapa faktor yang menjadi landasan dan membuat kami memilih saham NVIDA. Berikut adalah latar belakang yang melandasi pemilihan ini:
- Dominasi di Industri Chip AI
NVIDIA Corporation tidak hanya merupakan pemimpin, tetapi juga menggambarkan dominasi dalam industri chip kecerdasan buatan (AI). Sejak pendiriannya pada tahun 1993, perusahaan ini telah memainkan peran kunci dalam menghadirkan solusi grafis dan komputasi untuk berbagai industri, termasuk gaming, visualisasi profesional, pusat data, dan otomotif. - Inovasi Chip H200 sebagai Pemacu Utama
Pengumuman pengembangan chip terbaru, H200, menambahkan dimensi baru pada potensi inovatif NVIDIA dalam mendukung perkembangan AI. Dengan peningkatan signifikan pada memori dan kecepatan, H200 memberikan daya ungkit yang kuat bagi aplikasi AI generatif dan model bahasa besar, menjadikannya subjek yang menarik untuk diteliti. - Peran Penting dalam Layanan AI Generatif
NVIDIA tidak hanya mendominasi pasar chip AI, tetapi juga memainkan peran penting dalam layanan AI generatif, termasuk layanan ChatGPT yang populer. Sebagai penyedia teknologi di balik berbagai aplikasi AI canggih, perusahaan ini memiliki dampak yang signifikan pada perkembangan dan kemajuan kecerdasan buatan. - Dukungan dari Cloud Service Providers Terkemuka
Keterlibatan NVIDIA dengan penyedia Cloud Service terkemuka seperti Amazon Web Services, Google Cloud, Microsoft Azure, dan Oracle Cloud Infrastructure menunjukkan pengakuan industri terhadap kualitas dan relevansi produk-produknya. Ini menambah nilai pada data saham NVIDIA sebagai subjek penelitian. - Performa Saham yang Mengesankan
Kinerja saham NVIDIA yang mengesankan, terutama dengan kenaikan lebih dari tiga kali lipat dalam nilai selama tahun ini, menarik perhatian sebagai indikator potensial untuk menganalisis dampak inovasi teknologinya terhadap nilai pasar saham.
Dengan kombinasi faktor-faktor ini, pemilihan data saham NVIDIA bukan hanya memungkinkan analisis yang mendalam tentang hubungan antara inovasi teknologi dan pergerakan harga saham, tetapi juga memberikan wawasan yang berharga tentang bagaimana NVIDIA terus memainkan peran krusial dalam perkembangan industri kecerdasan buatan.
1.2 Dataset
Dataset yang saya gunakan merupakan koleksi data harga saham historis
periode 1 Januari 2018 hingga 10 November 2023 dari 6 saham yakni saham
saingan NVDA seperti AMD, ARM,
AVGO, TSM, dan INTC.
Dataset ini memilki data :
- Open: yakni Harga saham pada awal periode perdagangan tertentu. Ini adalah harga saham pertama pada hari perdagangan tersebut.
- High: Harga tertinggi yang saham capai selama periode perdagangan tersebut. Ini mencerminkan harga tertinggi yang pembeli bersedia bayar selama hari tersebut.
- Low: Harga terendah yang saham capai selama periode perdagangan tersebut. Ini mencerminkan harga terendah yang penjual bersedia terima selama hari tersebut.
- Close: Harga saham pada akhir periode perdagangan tertentu. Ini adalah harga saham terakhir pada hari perdagangan tersebut.
- Adj Close (Adjusted Close): Harga penutup yang telah disesuaikan untuk memperhitungkan perubahan seperti pembagian saham atau dividen. Ini adalah harga penutup yang paling relevan untuk analisis jangka panjang, karena mencerminkan harga saham yang sebenarnya setelah penyesuaian.
- Volume: Volume perdagangan saham selama periode tertentu. Ini mencerminkan jumlah saham yang diperdagangkan selama hari perdagangan tersebut.
Kami akan menggunakan peubah Adj Close (Adjusted Close),
Karena sesuai dengan penjelasan diatas, peubah Adj Close
adalah yang paling sesuai untuk dianalisis dibandingkan peubah
lainnya.
1.3 Tujuan
Tujuan dari praktikum ini adalah untuk menganalisis pola perkiraan pergerakan harga tujuh saham teknologi terkemuka dengan harapan dapat memberikan rekomendasi kepada pembaca mengenai saham mana yang sebaiknya dipertimbangkan untuk dibeli atau diinvestasikan secara signifikan di antara tujuh perusahaan teknologi besar tersebut.
1.4 Data Preparation
1.4.1 Import Data
install_load('rio')
raw.data <- import("https://raw.githubusercontent.com/Zen-Rofiqy/STA1341-MPDW/main/Data/New/%40AAANTI%20Stock%20Prices.csv")
1.4.2 Data Checking
Cek Tipe data.
str(raw.data)
## 'data.frame': 7464 obs. of 8 variables:
## $ Name : chr "AMD" "AMD" "AMD" "AMD" ...
## $ Date : IDate, format: "2018-01-02" "2018-01-03" ...
## $ Open : num 10.4 11.6 12.1 12.2 12 ...
## $ High : num 11 12.1 12.4 12.2 12.3 ...
## $ Low : num 10.3 11.4 12 11.7 11.8 ...
## $ Close : num 11 11.6 12.1 11.9 12.3 ...
## $ Adj Close: num 11 11.6 12.1 11.9 12.3 ...
## $ Volume : int 44146300 154066700 109503000 63808900 63346000 62560900 52561200 38354900 47149300 42686600 ...
Semua data Karakter, harus diubah.
Cek Data kosong.
sum(is.na(raw.data))
## [1] 0
Tidak ada data kosong.
1.4.3 Penyesuaian Tipe Data
Semua tipe data masih berupa character. Harus diubah menjadi tipe data yang sesuai.
install_load('dplyr')
data <- raw.data %>%
mutate(
Date = as.Date(raw.data[, 2], format = "%m/%d/%y"), #Mengubah menjadi Date
across(3:ncol(raw.data), as.numeric) #Mengubah menjadi Numerik
)
str(data)
## 'data.frame': 7464 obs. of 8 variables:
## $ Name : chr "AMD" "AMD" "AMD" "AMD" ...
## $ Date : Date, format: "2018-01-02" "2018-01-03" ...
## $ Open : num 10.4 11.6 12.1 12.2 12 ...
## $ High : num 11 12.1 12.4 12.2 12.3 ...
## $ Low : num 10.3 11.4 12 11.7 11.8 ...
## $ Close : num 11 11.6 12.1 11.9 12.3 ...
## $ Adj Close: num 11 11.6 12.1 11.9 12.3 ...
## $ Volume : num 4.41e+07 1.54e+08 1.10e+08 6.38e+07 6.33e+07 ...
1.4.4 Rechecking Data
Cek kembali data kosong.
cat('Banyaknya Data Kosong', sum(is.na(data)))
## Banyaknya Data Kosong 0
1.4.5 Cek Periode Data
data2 <- data
install_load("lubridate")
dates <- as.Date(data2$Date)
# Buat rentang waktu mulai dari tanggal pertama hingga tanggal terakhir dalam data
full_date_range <- seq(min(dates), max(dates), by = "days")
# Bandingkan rentang waktu dengan tanggal yang ada dalam data
missing_dates <- setdiff(full_date_range, dates)
# Jika 'missing_dates' kosong, maka semua tanggal sudah ada dalam data
if (length(missing_dates) == 0) {
cat("Semua tanggal ada dalam data.\n")
} else {
cat("Tanggal yang tidak ada dalam data sebanyak", length(missing_dates),
"\nAtau sebanyak", length(missing_dates) * 6,
"Data Hilang dari ke-6 perusahaan yang ada")
}
## Tanggal yang tidak ada dalam data sebanyak 667
## Atau sebanyak 4002 Data Hilang dari ke-6 perusahaan yang ada
Inputasi Data
install_load('purrr')
# Fungsi untuk mengisi data yang hilang
fill_missing_data <- function(name) {
data_filtered <- data2 %>%
filter(Name == name)
full_date_range <- seq(min(data2$Date), max(data2$Date), by = "days")
data_frame_template <- data.frame(Date = full_date_range)
# Menambahkan kolom "Name" sesuai dengan perusahaan yang diproses
data_frame_template$Name <- name
data_filled <- merge(data_frame_template, data_filtered,
by = c("Date", "Name"), all.x = TRUE)
return(data_filled)
}
# Menggunakan purrr::map untuk memproses setiap nama perusahaan
filled_data_list <- map(unique(data2$Name), fill_missing_data)
# Gabungkan data-data yang telah diisi menjadi satu data frame
final_data <- data.frame()
final_data <- do.call(rbind, filled_data_list)
# Urutkan data berdasarkan "Name" terlebih dahulu, kemudian "Date"
final_data <- final_data %>%
dplyr::select(1, 2, 7) %>%
arrange(Name, Date)
#Input Data Hilang
install_load('imputeTS')
data <- na_interpolation(final_data$`Adj Close`, option = "linear")
install_load('ggplot2')
chart <- ggplot_na_imputations(final_data$`Adj Close`, data)
chart
#Export Chart
ggsave("00_Imputasi Data.png", chart, path = export.chart,
dpi = 300, height = 12, width = 23)
data <- final_data %>% dplyr::select(Name, Date) %>%
mutate(`Adj Close` = data)
Data sudah di imputasi.
Cek ukuran data
cat("Ukuran data awal adalah", nrow(data2), "(", nrow(data2)/6, "Jika persaham)",
"\nBanyaknya peride data yang hilang", (nrow(data)-nrow(data2))/6,
"(PerSaham.", (nrow(data)-nrow(data2)) , "Jika semua)",
"\nUkuran data yang baru seharusnya =",
nrow(data2) + (nrow(data)-nrow(data2)) ,
"\nIni sudah sesuai dengan Ukuran data input yakni =", nrow(data),
"(", nrow(data)/6, "jika persaham)",
"\nJadi ada sekitar", (nrow(data)-nrow(data2))/nrow(data)*100,"% Data hilang")
## Ukuran data awal adalah 7464 ( 1244 Jika persaham)
## Banyaknya peride data yang hilang 906 (PerSaham. 5436 Jika semua)
## Ukuran data yang baru seharusnya = 12900
## Ini sudah sesuai dengan Ukuran data input yakni = 12900 ( 2150 jika persaham)
## Jadi ada sekitar 42.13953 % Data hilang
Data sudah benar, siap untuk dianalisis.
Jika dalam rentang tahun 2021 hingga 2023 akan seperti :
coba <- data %>%
filter(Date >= as.Date("2021-01-01") & Date <= as.Date("2023-11-10"))
coba2 <- data2 %>%
filter(Date >= as.Date("2021-01-01") & Date <= as.Date("2023-11-10"))
cat("Ukuran data awal adalah", nrow(coba2), "(", nrow(coba2)/6, "Jika persaham)",
"\nBanyaknya peride data yang hilang", (nrow(coba)-nrow(coba2))/6,
"(PerSaham.", (nrow(coba)-nrow(coba2)) , "Jika semua)",
"\nUkuran data yang baru seharusnya =",
nrow(coba2) + (nrow(coba)-nrow(coba2)) ,
"\nIni sudah sesuai dengan Ukuran data input yakni =", nrow(coba),
"(", nrow(coba)/6, "jika persaham)",
"\nJadi ada sekitar", (nrow(coba)-nrow(coba2))/nrow(coba)*100,"% Data hilang")
## Ukuran data awal adalah 3642 ( 607 Jika persaham)
## Banyaknya peride data yang hilang 437 (PerSaham. 2622 Jika semua)
## Ukuran data yang baru seharusnya = 6264
## Ini sudah sesuai dengan Ukuran data input yakni = 6264 ( 1044 jika persaham)
## Jadi ada sekitar 41.85824 % Data hilang
1.4.6 Data Cleaned
install_load('DT')
datatable(data, filter = 'top',
options = list(pageLength = 5))
2 Eksplorasi Data
col.amd <- c("#D02A49", "#F3AEAA"); col.arm <- c("#4EC2C1", "#B8E0DC")
col.avgo <- c("#DC4F26", "#FCCFBC"); col.nvda <- c("#36B450", "#ACD694")
col.tsm <- c("#EDB64E", "#F6DEB3"); col.intc <- c("#4B75BA", "#9BC6EA")
cols <- c("AMD"=col.amd[1], "test_amd"=col.amd[2],
"ARM"=col.arm[1], "test_arm"=col.arm[2],
"AVGO"=col.avgo[1], "test_avgo"=col.avgo[2],
"NVDA"=col.nvda[1], "test_mvda"=col.nvda[2],
"TSM"=col.tsm[1], "test_tsm"=col.tsm[2],
"INTC"=col.intc[1], "test_intc"=col.intc[2])
2.1 Plot Time Series
install_load('ggplot2','extrafont')
# font_import(); loadfonts() #Run ini sekali aja
theme.ts <- list(
theme(legend.position = "none",
axis.text.x = element_text(hjust = 1,
margin = margin(b = 10, t=20)),
axis.text.y = element_text(vjust = 0.5, face = "bold",
margin = margin(l = 20, r = 20)),
plot.title = element_text(hjust = 0.5, face = "bold"),
text = element_text(size = 30),
plot.subtitle = element_text(hjust = 0.5),
panel.background = element_rect(fill = 'transparent'),
plot.background = element_rect(fill='transparent', color=NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(linewidth = 1, colour = "black"))
)
theme.ts1 <- list(
theme(legend.position = "none",
axis.text.x = element_text(hjust = 1,
margin = margin(b = 10, t=20)),
axis.text.y = element_text(vjust = 0.5, face = "bold",
margin = margin(l = 50, r = 20)),
plot.title = element_text(hjust = 0.5, face = "bold"),
text = element_text(size = 30),
plot.subtitle = element_text(hjust = 0.5),
panel.background = element_rect(fill = 'transparent'),
plot.background = element_rect(fill='transparent', color=NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(linewidth = 1, colour = "black"))
)
2.1.1 AAANTI
Melihat keseluruhan Time Series data saham.
install_load('viridis','ggrepel')
#Plot
chart <-
ggplot(data, aes(x=Date, y=`Adj Close`, color=Name, alpha=Name)) + #Data
geom_line(aes(color=Name), linewidth=1.5) + #Timeseries
#Color
scale_color_manual(values = c(AMD=col.amd[1], ARM=col.arm[1],
AVGO=col.avgo[1], NVDA=col.nvda[1],
TSM=col.tsm[1], INTC=col.intc[1] )) +
scale_alpha_manual(values = c("AMD" = .25, "ARM" = .25, "AVGO" = .25,
"NVDA" = 1, "TSM" = .25, "INTC" = .25)) +
theme.ts + #THeme
labs(x = "\nPeriode (Tahun)", y='Harga Saham (USD)',
title = "Time Series NVDA",
subtitle = "Seperti apa sih pola deret waktu saham NVDA dan saingannya?\n") +
# Label / legend
geom_text_repel(
data=data[data$Date == max(data$Date),], #Posisi di ujung data
aes(color = Name, label = Name), #Warna garis & label saham
size = 8, #Ukuran text
nudge_x = 80, #Posisi Text (kanan 50)
hjust = 0, #Ujung
segment.size = 1, #Ukuran garis
segment.alpha = .75, #transparasi garis
segment.linetype = "dotted", #Time garis
box.padding = .4, #Biar label saham nggak dempetan
segment.curvature = -0.1, #biar garis mulus
segment.ncp = 8,
segment.angle = 60
) +
#Axis
coord_cartesian(clip = "off"
) +
scale_x_date( #Sumbu x
date_breaks = "1 year", # Menampilkan label setiap tahun
date_labels = "%Y", # Format label tahun
limits = c(as.Date(min(data$Date)),
as.Date(max(data$Date)) + 120)
#Tampilin lebih dari 20023-07-28 agar label saham bisa masuk
) +
scale_y_continuous( #Sumbu y
labels = scales::dollar_format(prefix = "$") #tambahin dolar
) +
annotate( #Buat nandain batas data
"text", x = as.Date(max(data$Date)), y = -35,
label = max(data$Date), size=6
) +
geom_vline( #Buat garis batas data
xintercept = as.numeric(as.Date( max(data$Date) )),
linetype = "dotted", color = "red")
chart
#Export Chart
ggsave("01_Time Series AAANTI.png", chart, path = export.chart,
dpi = 300, height = 12, width = 27)
Data saham berakhir pada tanggal 21 November 2023 dengan harga saham
NVDA merupakan saham tertinggi kedua. Jika dilihat dari
tahun 2019-2022, semua saham cenderung memiliki pola trend
naik. Lalu dari 2022-2023 polanya cenderung trend turun.
Dan 2023 keatas, polanya cenderung naik.
data2 <- data %>%
filter(Date >= as.Date("2021-01-01") & Date <= as.Date("2023-11-10"))
install_load('viridis','ggrepel')
#Plot
chart <-
ggplot(data2, aes(x=Date, y=`Adj Close`, color=Name, alpha=Name)) + #Data
geom_line(aes(color=Name), linewidth=1.5) + #Timeseries
#Color
scale_color_manual(values = c(AMD=col.amd[1], ARM=col.arm[1],
AVGO=col.avgo[1], NVDA=col.nvda[1],
TSM=col.tsm[1], INTC=col.intc[1] )) +
scale_alpha_manual(values = c("AMD" = .35, "ARM" = .35, "AVGO" = .35,
"NVDA" = 1, "TSM" = .35, "INTC" = .35)) +
theme.ts + #THeme
labs(x = "\nPeriode (Tahun)", y='Harga Saham (USD)',
title = "Time Series NVDA 2021-2023",
subtitle = "Seperti apa sih pola deret waktu saham NVDA dan saingannya?\n") +
# Label / legend
geom_text_repel(
data=data2[data2$Date == max(data2$Date),], #Posisi di ujung data
aes(color = Name, label = Name), #Warna garis & label saham
size = 8, #Ukuran text
nudge_x = 20, #Posisi Text (kanan 50)
hjust = 0, #Ujung
segment.size = 1, #Ukuran garis
segment.alpha = .75, #transparasi garis
segment.linetype = "dotted", #Time garis
box.padding = .4, #Biar label saham nggak dempetan
segment.curvature = -0.1, #biar garis mulus
segment.ncp = 8,
segment.angle = 60
) +
#Axis
coord_cartesian(clip = "off"
) +
scale_x_date( #Sumbu x
date_breaks = "1 year", # Menampilkan label setiap tahun
date_labels = "%Y", # Format label tahun
limits = c(as.Date(min(data2$Date)),
as.Date(max(data2$Date)) + 60)
#Tampilin lebih dari 20023-07-28 agar label saham bisa masuk
) +
scale_y_continuous( #Sumbu y
labels = scales::dollar_format(prefix = "$") #tambahin dolar
) +
annotate( #Buat nandain batas data
"text", x = as.Date(max(data2$Date)), y = -35,
label = max(data2$Date), size=6
) +
geom_vline( #Buat garis batas data
xintercept = as.numeric(as.Date( max(data2$Date) )),
linetype = "dotted", color = "red")
chart
#Export Chart
ggsave("01_Time Series AAANTI_2022-2023.png", chart, path = export.chart,
dpi = 300, height = 12, width = 27)
2.1.2 NVDA
nvda <- data2 %>%
filter(Name == "NVDA") # Filter data saham Amazon tahun 2022 ke atas
rownames(nvda) <- NULL
str(nvda)
## 'data.frame': 1044 obs. of 3 variables:
## $ Name : chr "NVDA" "NVDA" "NVDA" "NVDA" ...
## $ Date : Date, format: "2021-01-01" "2021-01-02" ...
## $ Adj Close: num 130 131 131 131 134 ...
datatable(nvda, filter = 'top',
options = list(pageLength = 5))
Mengubah Ajd Close Menjadi Time series.
nvda.ts <- ts(nvda[,3])
Ringkasan Data Ajd CLose.
summary(nvda.ts)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 112.2 158.1 206.4 237.1 278.5 493.5
install_load('ggtext')
min_value <- min(nvda$`Adj Close`)
min_date <- nvda$Date[which.min(nvda$`Adj Close`)]
min_perc <- (which.min(nvda$`Adj Close`) / nrow(nvda)) * 100
max_value <- max(nvda$`Adj Close`)
max_date <- nvda$Date[which.max(nvda$`Adj Close`)]
max_perc <- (which.max(nvda$`Adj Close`) / nrow(nvda)) * 100
chart <-
ggplot(nvda, aes(x=Date, y=`Adj Close`)) +
geom_line(aes(color=Name), linewidth=2) +
scale_color_manual(values = col.nvda[1]) +
labs(x = "\nPeriode (Tahun)", y='Saham Harga penutup',
title = "Time Series Saham NVIDIA",
subtitle = "Seperti apa sih pola deret waktu saham NVIDIA?\n") +
theme(legend.position = "none") +
theme.ts1 +
#Ekstras
#Titik terendah
geom_segment(aes(x = min_date,
xend = min_date,
y = min(nvda$`Adj Close`)* 1.01,
yend = max(nvda$`Adj Close`)*40/100),
arrow = arrow(type = "closed", length = unit(0.1, "inches")),
lineend = "round", color = "#D02A49", size=1.5) +
geom_richtext(
data = data.frame(x = min_date,
y = max(nvda$`Adj Close`)*40/100,
label = paste0("Titik Terendah") ),
aes(x, y, label = label), size = 7, color = "white",
fill = "#D02A49", box.color = "white", parse = TRUE
) +
geom_text(aes(x = min_date-1*35, y = max(`Adj Close`)*40/100, label =
paste0(min_date)),
vjust = -1.5, hjust = 0, size = 7, color = "grey30") +
#Titik tertinggi
geom_segment(aes(x = max_date,
xend = max_date,
y = max(nvda$`Adj Close`)* .99,
yend = max(nvda$`Adj Close`)*70/100),
arrow = arrow(type = "closed", length = unit(0.1, "inches")),
lineend = "round", color = "#4B75BA", size=1.5) +
geom_richtext(
data = data.frame(x = max_date,
y = max(nvda$`Adj Close`)*70/100,
label = paste0("Titik Tertinggi") ),
aes(x, y, label = label), size = 7, color = "white",
fill = "#4B75BA", box.color = "white", parse = TRUE
) +
geom_text(aes(x = max_date-1*35, y = max(`Adj Close`)*60/100, label =
paste0(max_date)),
vjust = -1.5, hjust = 0, size = 7, color = "grey30") +
annotate( #Buat nandain batas data
"text", x = as.Date(max(data2$Date)) -32, y = 50,
label = max(data2$Date), size=6
) +
geom_vline( #Buat garis batas data
xintercept = as.numeric(as.Date( max(data2$Date) )),
linetype = "dotted", color = "red", linewidth=1.5)
chart
#Export Chart
ggsave("02_Time Series NVIDIA.png", chart, path = export.chart,
dpi = 300, height = 12, width = 27)
Berdasarkan grafik deret waktu yang disajikan, terlihat bahwa dari
tahun 2021-2022 terjadi kenaikan harga saham/tren
naik. Sedangkan 2022-2023 terjadi penurunan harga
saham/tren turun. Dan pada tahun 2023 akhir tahun
terjadi kenaikan harga saham yang drastis/tren
naik.
2.1.3 Seasonal
install_load("forecast")
## Warning: package 'forecast' was built under R version 4.2.3
ggseasonplot(nvda.ts)
## Error in ggseasonplot(nvda.ts): Data are not seasonal
Terlihat bahwa data tidak memiliki pola musiman.
2.2 Data Train vs Test
#membagi 80% data latih (training) dan 20% data uji (testing)
train_nvda <- nvda[1: round(nrow(nvda) *78/100),]
test_nvda <- nvda[round(nrow(nvda) *78/100 +1): nrow(nvda),]
train_nvda.ts <- ts(train_nvda[,3])
test_nvda.ts <- ts(test_nvda[,3])
chart <-
ggplot() +
#Label Data Asli
annotate( "rect", alpha=.1, fill="#4EC2C1",
xmin=as.Date(min(data2$Date)),
xmax=as.Date(data2$Date[nrow(train_nvda)]),
ymin=max(nvda$`Adj Close`) * 1.25, ymax=Inf ) +
annotate( "text", color="#4EC2C1",
x = as.Date(data2$Date[ceiling(nrow(train_nvda)/2)]),
y = max(nvda$`Adj Close`) * 1.3,
label = "Data Latih", size=10) +
annotate( "text", color="#4EC2C1",
x = as.Date(data2$Date[ceiling(nrow(train_nvda)/2)]),
y = max(nvda$`Adj Close`) * 1.2,
label = paste0(nrow(train_nvda)," Hari (",
round(nrow(train_nvda)/nrow(nvda)*100, 1),
"%)"), size=10) +
#Label Data Ramal
annotate( "rect", alpha=.1, fill="violetred",
xmin=as.Date(data2$Date[nrow(train_nvda)]),
xmax=as.Date(max(data2$Date)),
ymin=max(nvda$`Adj Close`) * 1.25, ymax=Inf ) +
annotate( "text", color="violetred",
x = as.Date(data2$Date[ceiling( .89*length(data2$Date)/6)]) ,
y = max(nvda$`Adj Close`) * 1.3,
label = "Data Uji", size=10) +
annotate( "text", color="violetred",
x = as.Date(data2$Date[ceiling( .89*length(data2$Date)/6)]) ,
y = max(nvda$`Adj Close`) * 1.2,
label = paste0(nrow(test_nvda)," Hari (",
round(nrow(test_nvda)/nrow(nvda)*100, 1),
"%)"), size=10) +
geom_vline( #Buat garis batas data
xintercept = as.Date(data2$Date[nrow(train_nvda)]) ,
linetype = "dotted", color = "black", linewidth = 2) +
#NVDA
geom_line(data = train_nvda, linewidth=2,
aes(x = Date, y = `Adj Close`, col = "NVDA")) +
geom_line(data = test_nvda, linewidth=2,
aes(x = Date, y = `Adj Close`, col = "test_nvda")) +
geom_point(data = tail(train_nvda, 1), alpha = .5,
aes(x = Date, y = `Adj Close`), stroke=2,
size = 15, shape = 21, color = "black", fill="violetred") +
scale_colour_manual(values = c("NVDA"=col.nvda[1], "test_nvda"=col.nvda[2])) +
theme.ts + #THeme
labs(x = "\nPeriode (Tahun)", y='Harga Saham (USD)',
title = "Time Series Saham NVDA",
subtitle = "Pembagian Data Latih dan data Uji\n") +
#Axis
coord_cartesian(clip = "off"
) +
scale_y_continuous( #Sumbu y
labels = scales::dollar_format(prefix = "$") #tambahin dolar
) +
#Bagi Data Train & Test
geom_segment(aes(x = as.Date(max(train_nvda$Date)),
xend = as.Date(max(train_nvda$Date)),
y = max(nvda$`Adj Close`)* .58,
yend = max(nvda$`Adj Close`)*.7),
arrow = arrow(type = "closed", length = unit(0.1, "inches")),
lineend = "round", color = "#4B75BA", size=1.5) +
geom_richtext(
data = data.frame(x = as.Date(max(train_nvda$Date)),
y = max(nvda$`Adj Close`)*.7,
label = max(train_nvda$Date) ),
aes(x, y, label = label), size = 8, color = "white",
fill = "#4B75BA", box.color = "white", parse = TRUE
) +
geom_segment(aes(x = as.Date(max(data2$Date)),
xend = as.Date(max(data2$Date)),
y = max(nvda$`Adj Close`)* .95,
yend = max(nvda$`Adj Close`)*1.05),
arrow = arrow(type = "closed", length = unit(.1, "inches")),
lineend = "round", color = "#D02A49", size=1.5) +
geom_richtext(
data = data.frame(x = as.Date(max(data2$Date)),
y = max(nvda$`Adj Close`)*1.05,
label = max(data2$Date) ),
aes(x, y, label = label), size = 8, color = "white",
fill = "#D02A49", box.color = "white", parse = TRUE
)
chart
#Export Chart
ggsave("02_TS_NVDA_train-test.png", chart, path = export.chart,
dpi = 300, height = 12, width = 27)
Berdasarkan plot data deret waktu pada data latih (\(78\%\) dari data asli), terlihat bahwa data menunjukkan tren naik dan tren turun. Ini mengisyaratkan bahwa data latih tidak memenuhi kriteria stasioneritas dalam rataan maupun ragam. Di sisi lain, dalam plot data uji (\(22\%\) dari data asli), terlihat adanya tren yang melonjak naik dan kurangnya nilai tengah yang stabil. Ini juga menunjukkan bahwa data uji tidak stasioner dalam rataan.
3 Stasioneritas
3.1 Uji Stasioneritas
3.1.1 Plot ACF
install_load('tsibble','tseries')
acf(train_nvda.ts, main="", col=col.nvda[1], lwd=2)
mtext("NVDA", side=3, line=1, cex=2, font=2)
Berdasarkan plot ACF, terlihat bahwa plot ACF data train menurun secara perlahan (tails of slowly). Hal ini juga menjadi indikasi bahwa data tidak stasioner dalam rataan dan tidak membentuk gelombang sinus.
3.1.2 Uji ADF
tseries::adf.test(train_nvda.ts)
##
## Augmented Dickey-Fuller Test
##
## data: train_nvda.ts
## Dickey-Fuller = -1.6692, Lag order = 9, p-value = 0.7183
## alternative hypothesis: stationary
if (tseries::adf.test(train_nvda.ts)[["p.value"]] < 0.05) {
cat("Karena p-value < 0.05, Maka Rataan Stasioner")
} else {
cat("Karena p-value > 0.05, Maka Rataan Tidak Stasioner")
}
## Karena p-value > 0.05, Maka Rataan Tidak Stasioner
\(H_0\) : Data tidak stasioner dalam rataan
\(H_1\) : Data stasioner dalam rataan
Berdasarkan uji ADF tersebut, p-value lebih besar dari taraf nyata \(5\%\) sehingga tak tolak \(H_0\) dan menandakan bahwa data tidak stasioner dalam rataan. Hal ini sesuai dengan hasil eksplorasi menggunakan plot time series dan plot ACF, sehingga ketidakstasioneran model kedepannya harus ditangani.
3.1.3 Plot Box-Cox
install_load('MASS')
index <- seq(1:nrow(train_nvda)) #Beda titik potong
bc =boxcox(train_nvda.ts~index, lambda=seq(-2, 4, by=.01))
title("NVDA", cex.main = 2)
lambda_nvda <- bc$x[which.max(bc$y)] #Nilai Rounded Lambda
sk_nvda <- bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)] #SK
Gambar di atas menunjukkan bahwa dari rentang nilai \(\lambda\) -2 hingga 4, pada selang yang dibuatnya tidak memuat nilai \(\lambda=1\), sehingga ragam tidak stasioner.
3.1.3.1 Lambda dan Selang Kepercayaan
lsk <- data.frame(
Lambda = c(lambda_nvda ),
Batas.Bawah = c( min(sk_nvda) ),
Batas.Atas = c(max(sk_nvda) )
) %>%
mutate(across(where(is.numeric), ~round(., 2)))
lsk$Keterangan <- apply(lsk, 1, function(row) {
ifelse(row["Batas.Bawah"] <= 1 && row["Batas.Atas"] >= 1,
"Ragam Stasioner", "Ragam Tidak Stasioner")
})
rownames(lsk) <- c("NVDA")
datatable(lsk, filter = "none")
Hal ini dibuktikan dengan tabel nilai rounded value (\(\lambda\)) dan selang kepercayaan diatas
yang menyatakan bahwa memuat nilai satu di dalam selangnya. Ini
mengindikasikan bahwa NVDA memerlukan transformasi data
karena tidak memiliki sifat stasioner dalam
ragam.
3.2 Penanganan Ketidakstasioneran Data
3.2.1 Dalam Ragam: Transformasi
Mengutip dari laman r-coder.com, formula transformasi data disesuaikan dengan nilai rounded value (\(\lambda\)) optimumnya.
| \(\mathbf{\lambda}\) | Transformasi | Transformasi Balik |
|---|---|---|
| -2 | \(1/x^2\) | \(1/\sqrt{x}\) |
| -1 | \(1/x\) | \(1/x\) |
| -0.5 | \(1/\sqrt{x}\) | \(1/x^2\) |
| 0 | \(\log{(x)}\) | \(e^x\) |
| 0.5 | \(\sqrt{x}\) | \(x^2\) |
| 1 | \(x\) | \(x\) |
| 2 | \(x^2\) | \(\sqrt{x}\) |
Sehingga transformasi yang digunakan pada data saham NVDIA adalah sebagai berikut.
lamb <- c(-2, -1, -.5, 0, .5, 1, 2)
# Fungsi untuk mengonversi lambda ke transformasi
trans <- function(lambda) {
lambda <- lamb[sapply(lambda, function(y) which.min(abs(y - lamb))) ]
case_when(
lambda == -2 ~ "1/x^2",
lambda == -1 ~ "1/x",
lambda == -0.5 ~ "1/sqrt(x)",
lambda == 0 ~ "log(x)",
lambda == 1 ~ "x",
lambda == 0.5 ~ "sqrt(x)",
lambda == 2 ~ "x^2",
TRUE ~ NA_character_
)
}
# Fungsi untuk mengonversi lambda ke transformasi balik
trans_balik <- function(lambda) {
lambda <- lamb[sapply(lambda, function(y) which.min(abs(y - lamb))) ]
case_when(
lambda == -2 ~ "1/sqrt(x)",
lambda == -1 ~ "1/x",
lambda == -0.5 ~ "1/x^2",
lambda == 0 ~ "exp(x)",
lambda == 1 ~ "x",
lambda == 0.5 ~ "x^2",
lambda == 2 ~ "sqrt(x)",
TRUE ~ NA_character_
)
}
# Menambahkan kolom Transformasi dan Trans Balik
tran <- lsk %>%
mutate(
Transformasi = trans(Lambda),
`Transformasi Balik` = trans_balik(Lambda)
) %>% dplyr::select(Lambda, Transformasi, `Transformasi Balik`)
datatable(tran, filter = 'none')
3.2.1.1 Transformasi
#Transformasi
train_nvda.ts.new <- 1/sqrt(train_nvda.ts)
#Plot BoxCox
index <- seq(1:nrow(train_nvda))
bc <-boxcox(train_nvda.ts.new ~index, lambda=seq(-2, 4, by=.01))
title("NVDA", cex.main=2)
lambda_nvda.new <- bc$x[which.max(bc$y)] #Nilai Rounded Lambda
sk_nvda.new <- bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)] #SK
lsk.new <- data.frame(
Lambda = c(lambda_nvda.new ),
Batas.Bawah = c( min(sk_nvda.new) ),
Batas.Atas = c( max(sk_nvda.new) )
) %>%
mutate(across(where(is.numeric), ~round(., 2)))
lsk.new$Keterangan <- apply(lsk.new, 1, function(row) {
ifelse(row["Batas.Bawah"] <= 1 && row["Batas.Atas"] >= 1,
"Ragam Stasioner", "Ragam Tidak Stasioner")
})
rownames(lsk.new) <- c("NVDA")
datatable(lsk.new, filter = 'none')
Setelah melalui proses transformasi, gambar dan tabel di atas mengindikasikan bahwa setiap rangkaian data sekarang mencakup nilai satu dalam selangnya, yang menunjukkan bahwa data tersebut sekarang sudah stasioner dalam ragam.
3.2.1.2 Plot Transformasi
ragam <- data.frame(
Date = train_nvda %>% dplyr::select(Date) ,
`Asli` = c(as.vector(train_nvda.ts) ),
`Trans` = c(as.vector(train_nvda.ts.new) ) *3000,
Name = rep(c("NVDA"),
each = nrow(train_nvda) )
)
colnames(ragam) <- c("Date", "Asli", "Trans", "Name")
chart <-
ggplot() +
# Akhir data
geom_vline( #Buat garis batas data
xintercept = as.numeric(as.Date( max(train_nvda$Date) )),
linetype = "dotted", color = "red") +
#Time Series
#Data Asli
geom_line(data = ragam, linewidth=2,
aes(x = Date, y = `Asli`, col = "Asli")) +
#Data Transformasi
geom_line(data = ragam, linewidth=2,
aes(x = Date, y = `Trans`, col = "Trans")) +
scale_colour_manual(values = c("Asli"=col.nvda[2], "Trans"=col.nvda[1])) +
theme.ts + #THeme
labs(x = "\nPeriode (Tahun)", y='Harga Saham (USD)',
title = "Plot Transformasi NVDA",
subtitle = "Perbandingan plot Transformasi dengan Plot Asli\n") +
#Axis
coord_cartesian(clip = "off"
) +
scale_x_date( #Sumbu x
date_breaks = "1 year", # Menampilkan label setiap tahun
date_labels = "%Y", # Format label tahun
limits = c(as.Date(min(ragam$Date)),
as.Date(max(ragam$Date)) + 80)
) +
scale_y_continuous( #Sumbu y
labels = scales::dollar_format(prefix = "$") #tambahin dolar
) +
ylim(min(ragam$Asli), max(ragam$Asli)) +
#Data asli
geom_segment(aes(x = as.Date(max(train_nvda$Date)),
xend = as.Date(max(train_nvda$Date))+60,
y = train_nvda$`Adj Close`[nrow(train_nvda)],
yend = train_nvda$`Adj Close`[nrow(train_nvda)]),
arrow = arrow(type = "closed", length = unit(.1, "inches")),
lineend = "round", color = "gray60", size=1.5) +
geom_richtext(
data = data.frame(x = as.Date(max(train_nvda$Date))+60,
y = train_nvda$`Adj Close`[nrow(train_nvda)],
label = "Data Asli" ),
aes(x, y, label = label), size = 8, color = "white",
fill = "gray60", box.color = "white", parse = TRUE
) +
# Data Transformasi
geom_segment(aes(x = as.Date(max(train_nvda$Date)),
xend = as.Date(max(train_nvda$Date)) + 60,
y = ragam$Trans[nrow(ragam)],
yend = ragam$Trans[nrow(ragam)] ),
arrow = arrow(type = "closed", length = unit(.1, "inches")),
lineend = "round", color = "#D02A49", size=1.5) +
geom_richtext(
data = data.frame(x = as.Date(max(train_nvda$Date)) +60,
y = ragam$Trans[nrow(ragam)],
label = "Transformasi" ),
aes(x, y, label = label), size = 8, color = "white",
fill = "#D02A49", box.color = "white", parse = TRUE
) +
#Tanggal akhir
geom_segment(aes(x = as.Date(max(train_nvda$Date)),
xend = as.Date(max(train_nvda$Date)),
y = -Inf,
yend = max(train_nvda$`Adj Close`)*.15),
arrow = arrow(type = "closed", length = unit(.1, "inches")),
lineend = "round", color = "gray60", size=1.5) +
geom_richtext(
data = data.frame(x = as.Date(max(train_nvda$Date)),
y = max(train_nvda$`Adj Close`)*.15,
label = max(train_nvda$Date) ),
aes(x, y, label = label), size = 8, color = "white",
fill = "gray60", box.color = "white", parse = TRUE
)
chart
#Export Chart
ggsave("03_Stas_4_Transform.png", chart, path = export.chart,
dpi = 300, height = 12, width = 27)
3.2.2 Dalam Rataan: Differencing
Differencing (diferensiasi) adalah suatu proses statistik yang digunakan untuk mengatasi ketidakstasioneran dalam data rata-rata atau tren data. Tujuannya adalah untuk membuat data menjadi lebih stasioner dengan mengurangi atau menghilangkan tren atau pola waktu yang mungkin ada dalam data tersebut.
#Diff
train_nvda.diff <- diff(train_nvda.ts.new, differences = 1)
#Plot
plot.ts(train_nvda.diff, lty=1, xlab="Periode (Hari)",
col = col.nvda[1], lwd = 2,
main="Plot Differencing", cex.main = 2)
Berdasarkan plot data deret waktu yang sudah di Differencing, terlihat bahwa seluruh data sudah stasioner dalam rataan ditandai dengan data bergerak pada nilai tengah tertentu (tidak terdapat trend ataupun musiman pada data).
3.3 Uji Ulang
Akan dilakukan uji ulang untuk mengecek kestasioneran data dalam rataan.
3.3.1 Plot ACF
stats::acf(train_nvda.diff, main="", col=col.nvda[1], lwd=2, ylim=c(-.1,.1))
mtext("NVDA", side=3, line=1, cex=2, font=2)
Berdasarkan ketujuh plot ACF tersebut, dapat dilihat bahwa lag nya sudah tidak tails off slowly (trennya telah berakhir dan korelasi antara pengamatan pada lag-lag yang lebih jauh telah berkurang secara signifikan), sehingga data dinyatakan stasioner dalam rataan. Ini menunjukkan bahwa masalah ketidakstasioneran dalam data telah berhasil diatasi.
3.3.2 Uji ADF
tseries::adf.test(train_nvda.diff)
##
## Augmented Dickey-Fuller Test
##
## data: train_nvda.diff
## Dickey-Fuller = -8.3563, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary
if (tseries::adf.test(train_nvda.diff)[["p.value"]] < 0.05) {
cat("Karena p-value < 0.05, Maka Rataan Stasioner")
} else {
cat("Karena p-value > 0.05, Maka Rataan Tidak Stasioner")
}
## Karena p-value < 0.05, Maka Rataan Stasioner
\(H_0\) : Data tidak stasioner dalam rataan
\(H_1\) : Data stasioner dalam rataan
Berdasarkan uji ADF tersebut, p-value yang didapatkan yakni lebih kecil dari taraf nyata \(5\%\) sehingga tolak \(H_0\) atau data stasioner dalam rataan. Hal ini sesuai dengan hasil eksplorasi menggunakan plot time series dan plot ACF, sehingga dalam hal ini ketidakstasioneran data sudah berhasil ditangani dan dapat dilanjutkan ke pemodelan.
4 Model
4.1 Model Tentatif
Identifikasi model tentatif dapat dilakukan dengan melihat pada plot ACF, PACF dan Tabel EACF. Mengutip dari postingan pada laman medium.com Panduan dasar untuk menginterpretasikan plot ACF dan PACF adalah sebagai berikut :
- Cari pola tail off pada ACF atau PACF.
- Jika tail off terjadi pada ACF → model AR → cut off pada PACF akan memberikan nilai p untuk AR(p).
- Jika tail off terjadi pada PACF → model MA → cut off pada ACF akan memberikan nilai q untuk MA(q).
- Jika terjadi tail off pada kedua ACF dan PACF → model ARMA.
4.1.1 Plot ACF & PACF
par(mfrow=c(1,2))
stats::acf(train_nvda.diff, main="", col=col.nvda[1], lwd=2, ylim=c(-.1,.1))
mtext("NVDA", side=3, line=1, cex=2, font=2)
stats::pacf(train_nvda.diff, main="", col=col.nvda[1], lwd=2, ylim=c(-.1,.1))
mtext("NVDA", side=3, line=1, cex=2, font=2)
Referensi : Significance level of ACF and PACF in R
acf.pacf <- function(dt){ #Fungsi Deteksi Cut off plot ACF & PACF
N <- acf(dt, plot = F)[["n.used"]] + 1 #Banyaknya Data
#Garis Signifikan
sl <- c(-1,1)*(exp(2*1.96/sqrt(N-3))-1)/(exp(2*1.96/sqrt(N-3))+1)
acf <- data.frame( lag = acf(dt, plot=F)[["lag"]],
acf = acf(dt, plot=F)[["acf"]] )
cut.off1 <- acf %>%
filter(acf <= min(sl) | acf >= max(sl)) %>%
dplyr::select(lag); colnames(cut.off1) <- "acf"
pacf <- data.frame(lag = pacf(dt, plot=F)[["lag"]],
acf = pacf(dt, plot=F)[["acf"]] )
cut.off2 <- pacf %>%
filter(acf <= min(sl) | acf >= max(sl)) %>%
dplyr::select(lag); colnames(cut.off2) <- "pacf"
cut.off <- merge(cut.off1, cut.off2, by = 0, all = TRUE)[-1]
return(cut.off)
}
data_frames <- list(
train_nvda.diff
)
name <- "NVDA"
acf.pacf.tab <- data.frame()
for(i in 1:length(name)){
dex <- data_frames[[i]] %>% acf.pacf() %>%
rename(!!paste0(name[i], ".a") := acf,
!!paste0(name[i], ".p") := pacf)
acf.pacf.tab <- merge(acf.pacf.tab, dex, by=0, all=TRUE)[-1]
}
datatable(acf.pacf.tab, filter = 'none')
Berdasarkan plot tersebut, terlihat bahwa plot ACF cut off pada lag 0. Sedangkan pada plot PACF cut off pada lag 14. Sehingga model tidak dapat di identifikasi dengan plot ACF.
4.1.2 Tabel EACF
Referensi : [University of South Carolina] Chapter 6: Model Specification for Time Series
Mengutip dari referensi diatas, Identifikasi model menggunakan Tabel EACF (Extended Autocorrelation Function) untuk proses ARMA(p, q) seharusnya secara teoritis memiliki pola segitiga nol, dengan nol di sudut kiri atas muncul pada baris ke-p dan kolom ke-q (dengan label baris dan kolom dimulai dari 0). Pola segitiga ini merupakan pola segitiga pertama dari setiap baris.
Sebagai contoh :
Dalam hal ini model tentatif yang terbentuk adalah ARIMA(0,1,1), ARIMA(1,1,1), ARIMA(2,1,2), ARIMA(3,1,3), dst.. Namun karena kemungkinannya sangat banyak, maka akan digunakan function agar menyingkat proses.
Pemilihan model ARIMA terbaik dilakukan ketika memiliki nilai AIC yang terkecil, Semua Parameter Signifikan, dan Tidak ada Parameter NA. Kami akan mencoba untuk menggunakan function pemodelam model ARIMA berdasarkan tabel EACF. Namun dengan melakukan 2 pemodelan perbaris AR(p) nya, agar sekalian overfitting.
install_load('TSA', 'forecast')
# Fungsi untuk menghitung model ARIMA dan menganalisis parameter
model_tentatif <- function(data, p_max, d, q_max, alpha=0.05, drift=FALSE) {
best_model <- NULL
best_aic <- Inf
eacf_result <- eacf(data) #Untuk Ektraksi & identifikasi model
models <- data.frame(Model = character(0),
AIC = numeric(0),
Signif = character(0),
Keterangan = character(0))
for (p in 0:p_max) {
overfit <- 0
for (q in 1:q_max) {
#Pola Matriks segitiga bawah
if (!is.na(eacf_result$symbol[p + 1, q ]) &&
!is.na(eacf_result$symbol[p + 1, q + 1]) &&
!is.na(eacf_result$symbol[p + 2, q + 1])) {
if (eacf_result$symbol[p + 1, q ] == "o" &&
eacf_result$symbol[p + 1, q + 1] == "o" &&
eacf_result$symbol[p + 2, q + 1] == "o") {
model <- Arima(data, order=c(p,d,q), method="ML",
include.drift=drift)
aic <- AIC(model)
# Mendapatkan nilai coef dari model
coeftest_result <- lmtest::coeftest(model)
# jika lebih kecil dari alpha, maka signifikan
significant_params <-
rownames(coeftest_result)[coeftest_result[, "Pr(>|z|)"] < alpha]
# jika lebih besar dari alpha, maka tidak signifikan
non_significant_params <-
rownames(coeftest_result)[coeftest_result[, "Pr(>|z|)"] > alpha]
# Keterangan signifikansi
if (length(significant_params) == 0) {
keterangan <- "Semua parameter tidak signifikan"
} else if (length(significant_params) == nrow(coeftest_result)) {
keterangan <- "Semua parameter signifikan"
} else {
keterangan <- paste("Parameter yang tidak signifikan adalah",
paste(non_significant_params, collapse = ", "))
}
models <- rbind(models,
data.frame(Model = paste("ARIMA(", p, ",", d, ",", q, ")",
sep = ""),
AIC = aic,
Signif = paste(significant_params,
collapse = ", "),
Keterangan = keterangan))
#Identifikasi Best Model
if (keterangan == "Semua parameter signifikan" &&
!any(is.na(significant_params))) {
if (aic < best_aic) {
best_model <- model
best_aic <- aic
}
}
#jika model ditemukan pada baris p. Maka lanjut ke p+1
if(overfit >= 1) break
overfit = overfit + 1
}
}
}
}
cat("\nModel ARIMA dengan AIC terkecil:\n")
print(best_model)
return(models)
}
model.tentaif_nvda <-
model_tentatif(train_nvda.diff, p_max = 6, d = 1, q_max = 12)
## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o o o o o o o o o o o o
## 1 x o o o o o o o o o o o o o
## 2 x x o o o o o o o o o o o o
## 3 x o x o o o o o o o o o o o
## 4 o x x x o o o o o o o o o o
## 5 x x x o x o o o o o o o o o
## 6 x x x x x x o o o o o o o o
## 7 x x o x x x o o o o o o o o
##
## Model ARIMA dengan AIC terkecil:
## Series: data
## ARIMA(0,1,1)
##
## Coefficients:
## ma1
## -0.9933
## s.e. 0.0078
##
## sigma^2 = 9.868e-07: log likelihood = 4460.68
## AIC=-8917.36 AICc=-8917.34 BIC=-8907.96
datatable(model.tentaif_nvda, filter = 'top',
options = list(pageLength = 7))
Berdasarkan pendugaan parameter di atas, nilai AIC terkecil \(-8917.35\) dimiliki oleh model
ARIMA(0,1,1) dan seluruh parameternya signifikan sehingga
model yang dipilih adalah model ARIMA(0,1,1).
model_nvda.da <- Arima(train_nvda.diff, order=c(0,1,1), method="ML")
4.2 Model ARIMA dengan Drift
Tren menentukan bagaimana perubahan deret waktu dari waktu ke waktu. Dalam model ARIMA, ini dapat dimodelkan menggunakan persamaan dasar berikut “a” disebut “intercept” dan “b” disebut “drift”. “Drift” adalah hanyalah kemiringan garis lurus. Selanjutnya akan dilakukan kombinasi model ARIMA dengan Intercept dan Drift.
4.2.1 1. Model Tanpa Intercept dan Tanpa Drift
Simpelnya ini merupakan model data train yang diberi differencing. Yakni memodelkan data train dengan ordo I dari ARIMA(p, d, q) sama dengan 1, sehingga kan membentuk model ARIMA(p, 1, q).
drift1_nvda <-
model_tentatif(train_nvda.ts, p_max = 6, d = 1, q_max = 12)
## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x x x x x x x x x x x x
## 1 o o o o o o o o o o o o o o
## 2 x o o o o o o o o o o o o o
## 3 x x o o o o o o o o o o o o
## 4 x x x o o o o o o o o o o o
## 5 x x x o o o o o o o o o o o
## 6 x x o x x o o o o o o o o o
## 7 x x x x x x o o o o o o o o
##
## Model ARIMA dengan AIC terkecil:
## Series: data
## ARIMA(4,1,4)
##
## Coefficients:
## ar1 ar2 ar3 ar4 ma1 ma2 ma3 ma4
## 0.4960 -0.7914 0.5091 -0.9828 -0.4804 0.7803 -0.5046 0.9699
## s.e. 0.0144 0.0120 0.0094 0.0172 0.0236 0.0193 0.0172 0.0240
##
## sigma^2 = 29.16: log likelihood = -2522.64
## AIC=5063.29 AICc=5063.51 BIC=5105.59
datatable(drift1_nvda, filter = 'top',
options = list(pageLength = 7))
Berdasarkan pendugaan parameter di atas, nilai AIC terkecil \(5063.53\) dimiliki oleh model
ARIMA(4,1,4) dan seluruh parameternya signifikan sehingga
model yang dipilih adalah model ARIMA(4,1,4).
model_nvda.drift1 <- Arima(train_nvda.ts, order=c(4,1,4), method="ML")
4.2.2 2. Model Dengan Intercept Tanpa Drift
Simpelnya ini merupakan model data differencing tanpa diberi differencing. Yakni memodelkan data differencing dengan ordo I dari ARIMA(p, d, q) sama dengan 0, sehingga kan membentuk model ARIMA(p, 0, q).
drift2_nvda <-
model_tentatif(train_nvda.diff, p_max = 6, d = 0, q_max = 12)
## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o o o o o o o o o o o o
## 1 x o o o o o o o o o o o o o
## 2 x x o o o o o o o o o o o o
## 3 x o x o o o o o o o o o o o
## 4 o x x x o o o o o o o o o o
## 5 x x x o x o o o o o o o o o
## 6 x x x x x x o o o o o o o o
## 7 x x o x x x o o o o o o o o
##
## Model ARIMA dengan AIC terkecil:
## NULL
datatable(drift2_nvda, filter = 'top',
options = list(pageLength = 7))
Kategori Model terbaik tidak ada yang terpenuhi,
yakni tidak ada nilai AIC terkecil dengan semua parameter nya signifikan
dan tidak ada NA di dalam nya. Sehingga Model Dengan Intercept Tanpa
Drift pada NVDA tidak ada.
model_nvda.drift2 <- NA
4.2.3 3. Model Tanpa Intercept Dengan Drift
Simpelnya ini merupakan model data train yang diberi differencing dan drift. Yakni memodelkan data train dengan ordo I dari ARIMA(p, d, q) sama dengan 1, sehingga kan membentuk model ARIMA(p, 1, q) + include.drift.
drift3_nvda <-
model_tentatif(train_nvda.ts, p_max = 6, d = 1, q_max = 12, drift=TRUE)
## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x x x x x x x x x x x x
## 1 o o o o o o o o o o o o o o
## 2 x o o o o o o o o o o o o o
## 3 x x o o o o o o o o o o o o
## 4 x x x o o o o o o o o o o o
## 5 x x x o o o o o o o o o o o
## 6 x x o x x o o o o o o o o o
## 7 x x x x x x o o o o o o o o
##
## Model ARIMA dengan AIC terkecil:
## NULL
datatable(drift3_nvda, filter = 'top',
options = list(pageLength = 7))
Kategori Model terbaik tidak ada yang terpenuhi,
yakni tidak ada nilai AIC terkecil dengan semua parameter nya signifikan
dan tidak ada NA di dalam nya. Sehingga Model Tanpa Intercept Dengan
Drift pada NVDA tidak ada.
model_nvda.drift3 <- NA
4.2.4 4. Model dengan Intercept dan dengan Drift
Simpelnya ini merupakan model data differencing tanpa diberi differencing namun diberi drift. Yakni memodelkan data differencing dengan ordo I dari ARIMA(p, d, q) sama dengan 0, sehingga kan membentuk model ARIMA(p, 0, q) + include.drift.
drift4_nvda <-
model_tentatif(train_nvda.diff, p_max = 6, d = 1, q_max = 12, drift=TRUE)
## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o o o o o o o o o o o o
## 1 x o o o o o o o o o o o o o
## 2 x x o o o o o o o o o o o o
## 3 x o x o o o o o o o o o o o
## 4 o x x x o o o o o o o o o o
## 5 x x x o x o o o o o o o o o
## 6 x x x x x x o o o o o o o o
## 7 x x o x x x o o o o o o o o
##
## Model ARIMA dengan AIC terkecil:
## NULL
datatable(drift4_nvda, filter = 'top',
options = list(pageLength = 7))
Kategori Model terbaik tidak ada yang terpenuhi,
yakni tidak ada nilai AIC terkecil dengan semua parameter nya signifikan
dan tidak ada NA di dalam nya. Sehingga Model dengan Intercept dan
dengan Drift pada NVDA tidak ada.
model_nvda.drift4 <- NA
4.3 Model Terbaik
model_terbaik <- function(model, rowname, mods) {
if (!all(is.na(model))) {
model_df <- data.frame(
Model = paste("ARIMA(",
paste(as.character(model[["call"]][["order"]][-1]),
collapse = ","), ")", sep = ""),
AIC = model[["aicc"]],
row.names = rowname
)
model_df$Keterangan <-
mods[which(mods$Model ==
model_df["Model"][rowname,]), "Keterangan"]
return(model_df)
} else {
return(data.frame(
Model = NA_character_,
AIC = NA_real_,
Keterangan = NA_character_,
row.names = rowname
))
}
}
best_model <- rbind(
model_terbaik(model_nvda.da, 'NVDA', model.tentaif_nvda),
model_terbaik(model_nvda.drift1, 'NVDA', drift1_nvda),
model_terbaik(model_nvda.drift2, 'NVDA', drift2_nvda),
model_terbaik(model_nvda.drift3, 'NVDA', drift3_nvda),
model_terbaik(model_nvda.drift4, 'NVDA', drift4_nvda)
) %>%
mutate(across(where(is.numeric), ~round(., 2)))
datatable(best_model, filter = 'top',
options = list(pageLength = 7))
Terlihat bahwa model paling terbaik adalah model tentatif tanpa
kombinasi drift. Sangat disayangkan, dari 4 kombinasi drift yang ada,
hanya kombinasi 1 yang memiliki kategori model terbaik. Sehingga model
yang akan digunakan adalah model_nvda.da.
4.4 Overfitting
Overfitting dilakukan dilakukan dengan menaikkan orde
AR(p) dan MA(q) dari model terpilih
ARIMA(0,1,1) untuk melihat apakah terdapat model lain yang
lebih baik dari model saat ini. Kandidat model overfitting adalah
ARIMA(1,1,1) dan ARIMA(0,1,2).
Perlu diingat bahwa sebelumnya saya menjelaskan bahwa fungsi
model_tentatif yang kami buat itu sudah sekaligus
overfitting, dan itu juga sudah membandingkan mana model terbaiknya.
Sehingga sebenarnya tidak perlu dilakukan overfitting. Namun saya akan
tunjukkan kembali.
datatable(model.tentaif_nvda, filter = 'top',
options = list(pageLength = 2))
Terlihat bahwa hanya terdapat model ARIMA(0,1,2) saja.
Itu dikarenakan model ARIMA(1,1,1) tidak tersedia di Tabel
EACF, sehingga perlu dilakukan manual.
Arima(train_nvda.diff, order=c(1,1,1), method="ML") %>% summary()
## Series: train_nvda.diff
## ARIMA(1,1,1)
##
## Coefficients:
## ar1 ma1
## 0.0040 -0.9934
## s.e. 0.0357 0.0081
##
## sigma^2 = 9.879e-07: log likelihood = 4460.68
## AIC=-8915.37 AICc=-8915.34 BIC=-8901.27
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -6.096225e-06 0.0009921142 0.0006814358 Inf Inf 0.7812342
## ACF1
## Training set -0.001297599
memiliki AIC yang lebih besar dan tidak semua
parameternya signifikan. Sehingga tetap model terbaiknya adalah
ARIMA(0,1,1).
arima.best_model_nvda <- model_nvda.da
arma.order <- c(
arima.best_model_nvda[["call"]][["order"]][[2]],
arima.best_model_nvda[["call"]][["order"]][[4]]
)
5 Analisis Sisaan
Model terbaik hasil identifikasi kemudian dicek asumsi sisaannya. Sisaan model ARIMA harus memenuhi asumsi normalitas, kebebasan sisaan, dan kehomogenan ragam. Diagnostik model dilakukan secara eksplorasi dan uji formal.
5.1 Eksplorasi Sisaan
#Eksplorasi
sisaan_nvda.da <- model_nvda.da$residuals
par(mfrow=c(2,2))
# QQ-Plot
qqnorm(sisaan_nvda.da)
qqline(sisaan_nvda.da, col = "dodgerblue3", lwd = 2.5)
# Plot sisaan
plot(c(1:length(sisaan_nvda.da)), sisaan_nvda.da,
main = "Plot Sisaan Model ARIMA",
xlab = "Periode",
ylab = "Nilai Sisaan")
abline(h = 0, col = "firebrick2", lty = 2, lwd=2.5)
#Plot ACF dan PACF
stats::acf(sisaan_nvda.da, ylim=c(-.1,.1))
stats::pacf(sisaan_nvda.da, ylim=c(-.1,.1))
Berdasarkan plot kuantil-kuantil normal, secara eksplorasi
ditunjukkan sisaan menyebar tidak normal ditandai
dengan tititk-titiknya cenderung tidak mengikuti garis \(45^{\circ}\). Kemudian dapat dilihat juga
lebar pita sisaan yang cenderung sama menandakan bahwa sisaan memiliki
ragam yang homogen. Plot ACF dan PACF sisaan
ARIMA(0,1,1) juga tidak signifikan pada 20 lag awal yang
menandakan saling bebas. Kondisi ini akan diuji lebih
lanjut dengan uji formal.
5.2 Uji Formal
formal <- data.frame(
Menguji = c("Kenormalan Siaan", "Kebebasan Sisaan",
"Homogenitas Sisaan", "Nilai Tengah = 0"),
Uji = c("Kolmogorov-Smirnov", "Ljung-Box",
"Ljung-Box Sisaan^2", "Uji-T"),
p.value = c(ks.test(sisaan_nvda.da, "pnorm")[["p.value"]],
Box.test(sisaan_nvda.da, type = "Ljung")[["p.value"]],
Box.test((sisaan_nvda.da)^2, type = "Ljung")[["p.value"]],
t.test(sisaan_nvda.da, mu=0, conf.level=.95)[["p.value"]])
)
formal$Keputusan[1] <- ifelse(formal$p.value[1] < .05, "Sisaan Tidak Menyebar Normal",
"Sisaan Menyebar Normal")
formal$Keputusan[2] <- ifelse(formal$p.value[2] < .05, "Sisaan Tidak Saling Bebas",
"Sisaan Saling Bebas")
formal$Keputusan[3] <- ifelse(formal$p.value[3] < .05, "Ragam Sisaan Tidak Homogen",
"Ragam Sisaan Homogen")
formal$Keputusan[4] <- ifelse(formal$p.value[4] < .05, "Nilai Tengah Sisaan \u2260 0",
"Nilai Tengah Sisaan = 0")
formal$p.value <- ifelse(round(formal$p.value, 2) == 0,
"< 2.2e-16",
round(formal$p.value, 2) )
datatable(formal, options = list(pageLength = 4))
5.2.1 Sisaan Menyebar Normal
Selain dengan eksplorasi, asumsi tersebut dapat diuji menggunakan uji formal. Pada tahapan ini uji formal yang digunakan untuk normalitas adalah uji Kolmogorov-Smirnov (KS). Hipotesis pada uji KS adalah sebagai berikut.
\(H_0\) : Sisaan menyebar normal
\(H_1\) : Sisaan tidak menyebar normal
ks.test(sisaan_nvda.da, "pnorm")
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: sisaan_nvda.da
## D = 0.49834, p-value < 2.2e-16
## alternative hypothesis: two-sided
if (ks.test(sisaan_nvda.da, "pnorm")[["p.value"]] < 0.05) {
cat("Karena p-value < 0.05, Maka Sisaan Tidak Menyebar Normal")
} else {
cat("Karena p-value > 0.05, Maka Sisaan Menyebar Normal")
}
## Karena p-value < 0.05, Maka Sisaan Tidak Menyebar Normal
Berdasarkan uji KS tersebut, didapat p-value yang sangat kecil yakni kurang dari \(2.2\times10^{-16}\). p-value nya kurang dari taraf nyata \(5\%\) sehingga tolak \(H_0\) dan menandakan bahwa sisaan tidak menyebar normal. Sesuai dengan qqplot sisaan.
5.2.2 Sisaan saling bebas
Selanjutnya akan dilakukan uji formal untuk kebebasan sisaan menggunakan uji Ljung-Box. Hipotesis yang digunakan adalah sebagai berikut.
\(H_0\) : Sisaan saling bebas
\(H_1\) : Sisaan tidak tidak saling bebas
Box.test(sisaan_nvda.da, type = "Ljung")
##
## Box-Ljung test
##
## data: sisaan_nvda.da
## X-squared = 0.005934, df = 1, p-value = 0.9386
if (Box.test(sisaan_nvda.da, type = "Ljung")[["p.value"]] < 0.05) {
cat("Karena p-value < 0.05, Maka Sisaan Tidak Saling Bebas")
} else {
cat("Karena p-value > 0.05, Maka Sisaan Saling Bebas")
}
## Karena p-value > 0.05, Maka Sisaan Saling Bebas
Berdasarkan uji Ljung-Box tersebut, p-value nya lebih besar dari taraf nyata \(5\%\) sehingga tak tolak \(H_0\) dan menandakan bahwa sisaan saling bebas. Yang berarti tidak ada autokorelasi. Hal ini sesuai dengan hasil eksplorasi menggunakan plot ACF dan PACF.
5.2.3 Sisaan homogen
Hipotesis yang digunakan untuk uji kehomogenan ragam adalah sebagai berikut.
\(H_0\) : Ragam sisaan homogen
\(H_1\) : Ragam sisaan tidak homogen
Box.test((sisaan_nvda.da)^2, type = "Ljung")
##
## Box-Ljung test
##
## data: (sisaan_nvda.da)^2
## X-squared = 4.2891, df = 1, p-value = 0.03836
if (Box.test((sisaan_nvda.da)^2, type = "Ljung")[["p.value"]] < 0.05) {
cat("Karena p-value < 0.05, Maka Sisaan Tidak Homogen")
} else {
cat("Karena p-value > 0.05, Maka Sisaan Homogen")
}
## Karena p-value < 0.05, Maka Sisaan Tidak Homogen
Berdasarkan uji Ljung-Box terhadap sisaan kuadrat tersebut, p-value nya lebih kecil dari taraf nyata \(5\%\) sehingga tolak \(H_0\) dan menandakan bahwa ragam sisaan tidak homogen.
5.2.4 Nilai tengah sama dengan nol
Terakhir, dengan uji-t, akan dicek apakah nilai tengah sisaan sama dengan nol. Hipotesis yang diujikan sebagai berikut.
\(H_0\) : nilai tengah sisaan sama dengan 0
\(H_1\) : nilai tengah sisaan tidak sama dengan 0
t.test(sisaan_nvda.da, mu=0, conf.level=.95)
##
## One Sample t-test
##
## data: sisaan_nvda.da
## t = -0.18269, df = 812, p-value = 0.8551
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -7.470093e-05 6.198009e-05
## sample estimates:
## mean of x
## -6.360421e-06
if (t.test(sisaan_nvda.da, mu=0, conf.level=.95)[["p.value"]] < 0.05) {
cat("Karena p-value < 0.05, Maka Nilai Tengah \u2260 0")
} else {
cat("Karena p-value > 0.05, Maka Nilai Tengah = 0")
}
## Karena p-value > 0.05, Maka Nilai Tengah = 0
Berdasarkan uji-t tersebut, p-value lebih besar dari taraf nyata \(5\%\) sehingga tak tolak \(H_0\) dan menandakan bahwa nilai tengah sisaan sama dengan nol.
6 Model ARCH-GARCH
Referensi : Measuring Volatility Using the ARCH Model, Volatility Modeling with R : ARCH and GARCH Models
6.1 ARCH
ARCH merupakan singkatan dari Auto-Regressive Conditional Teteroskedasticity.
- AutoRegressive (AR): Ini merujuk pada ide bahwa variansi (volatilitas) suatu variabel bergantung pada nilai-nilai sebelumnya dalam deret waktu. Dengan kata lain, variabilitas saat ini dipengaruhi oleh nilai-nilai sebelumnya.
- Conditional (Conditional): Ini menunjukkan bahwa perhitungan volatilitas dilakukan secara kondisional terhadap informasi sebelumnya dalam deret waktu. Dengan kata lain, variabilitas diukur dengan mempertimbangkan informasi kondisional terkini.
- Heteroskedasticity: Ini merujuk pada sifat ketidakseragaman (Heterogen/Tidak Homogen) dalam variabilitas. Dalam konteks ARCH, heteroskedasticity menunjukkan bahwa variabilitas suatu variabel tidak konstan, tetapi bervariasi sepanjang waktu berdasarkan informasi kondisional.
Jadi, AutoRegressive Conditional Heteroskedasticity (ARCH) adalah model statistik yang digunakan untuk mengukur dan memodelkan volatilitas dalam deret waktu finansial atau ekonometrik. Model ini mengasumsikan bahwa variabilitas harga atau hasil suatu variabel dapat bervariasi seiring waktu dan dipengaruhi oleh nilai-nilai sebelumnya dalam deret waktu.
Volatilitas sendiri merujuk pada tingkat fluktuasi atau variasi harga atau hasil suatu variabel dalam periode waktu tertentu. Volatilitas tinggi menunjukkan fluktuasi yang besar dalam nilai-nilai variabel tersebut, sementara volatilitas rendah menunjukkan fluktuasi yang lebih terbatas. Dalam konteks finansial, volatilitas sering kali dianggap sebagai ukuran risiko, karena fluktuasi harga yang besar dapat menunjukkan ketidakpastian atau risiko yang tinggi.
Model ARCH adalah model yang populer karena spesifikasi variansnya
dapat menangkap fitur umum dari deret waktu variabel keuangan. Model ini
berguna untuk memodelkan volatilitas dan khususnya
perubahan volatilitas seiring waktu. Mari kita lihat kembali pada plot
NVDA untuk mengevaluasi volatilitasnya.
#Plot
plot.ts(train_nvda.diff, lty=1, xlab="Periode (Hari)",
col = col.nvda[1], lwd = 2,
main="NVDA", cex.main = 2)
Seperti yang dapat Anda lihat, nilai dari deret ini berubah dengan cepat dari periode ke periode dalam suatu cara yang tampaknya tidak dapat diprediksi, yang dapat kita sebut sebagai volatil. Selain itu, ada periode-periode ketika perubahan besar diikuti oleh perubahan besar berikutnya dan periode di mana perubahan kecil diikuti oleh perubahan kecil berikutnya, merupakan bukti lain dari volatilitas.
6.1.1 Uji Efek ARCH
Uji Langrange multiplier (LM) sering digunakan untuk menguji keberadaan efek ARCH.
Untuk melakukan uji efek ARCH, kita harus
Mengestimasi persamaan rata-rata yang, dalam contoh ini, adalah \(r_t=\beta_0 + e_t\), dimana \(r\) =
train_nvda.diffMendapatkan residual yang diestimasi \(e^t\)
Mengestimasi \(\hat{e}^2_t = \gamma_0 + \gamma_1 \hat{e}^2_{t-1} + vt\)
dimana \(vt\) adalah suatu error acak.
Hipotesis
\(H_0\) : \(\gamma_1=0\) (Tidak ada efek ARCH)
\(H_1\) : \(\gamma_1 \ne0\) (Ada efek ARCH)
install_load('dynlm')
# Langkah 1: Mengestimasi persamaan rata-rata r = beta + error
byd.mean <- dynlm(train_nvda.diff ~ 1)
# Langkah 2: Mengambil residual dari model sebelumnya dan mengkuadratkannya
ehatsq <- ts(resid(byd.mean)^2)
# Langkah 3: Regresikan residual kuadrat pada residual kuadrat yang tertinggal satu lag
byd.arch <- dynlm(ehatsq ~ L(ehatsq), data = ehatsq)
summary(byd.arch)
##
## Time series regression with "ts" data:
## Start = 2, End = 813
##
## Call:
## dynlm(formula = ehatsq ~ L(ehatsq), data = ehatsq)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.235e-06 -9.005e-07 -7.837e-07 -2.000e-07 2.875e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.114e-07 8.638e-08 10.552 <2e-16 ***
## L(ehatsq) 7.268e-02 3.504e-02 2.074 0.0384 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.257e-06 on 810 degrees of freedom
## Multiple R-squared: 0.005282, Adjusted R-squared: 0.004054
## F-statistic: 4.301 on 1 and 810 DF, p-value: 0.0384
Kita dapat memeriksa efek ARCH dengan menggunakan fungsi
ArchTest() dari paket FinTS. Kita akan
menggunakan tingkat signifikansi \(\alpha=0.05\) untuk uji hipotesis nol kita.
Perlu diperhatikan bahwa ArchTest() memiliki default
lags=12.
install_load('FinTS')
lag <- 12
arch.effect <- data.frame(p.value = numeric(0),
Keterangan = character(0))
for(i in 1:lag){
p.value <- ArchTest(sisaan_nvda.da, lags=i, demean=TRUE)[["p.value"]]
arch.effect <- rbind(arch.effect,
data.frame(p.value = ifelse(round(p.value, 2) == 0,
sprintf("%.2e", p.value),
round(p.value, 2) ),
Keterangan = ifelse(p.value < 0.05,
"Signifikan", "Tidak signifikan")
) )
rownames(arch.effect)[i] <- paste0("Lag ", i)
}
datatable(arch.effect, filter = 'top',
options = list(pageLength = 6))
Dari Lag 1 sampai Lag 12 yang ditest, Terlihat bahwa Lag 2 memiliki p-value sebesar \(0.08\) yang lebih besar dari taraf nyata \(5\%\), sehingga tolak \(H_0\). Artinya Lag 2 Tidak sinifikan atau tidak ada efek ARCH(2). Selain Lag 2, Lag nya signifikan. Karena pada saat Lag 2 saja sudah tidak signifikan maka tidak perlu dilakukan model GARCH.
6.1.2 Model ARCH
install_load('rugarch')
create_spec <- function(i, j, arma_order, mean=mean) {
return(ugarchspec(
variance.model = list(model = "sGARCH", garchOrder = c(i, j)),
mean.model = list(armaOrder = arma_order,
#model tidak harus memasukkan komponen rata-rata (intercept)
include.mean = mean)
))
}
arch.garch <- function(data, pp, qq = 0, arma.order = c(0, 0), mean=FALSE) {
arch.model <- data.frame(Model = character(0),
Signif = numeric(0),
AIC = numeric(0),
Sig.par = character(0),
Non_Sig.par = character(0))
best_model <- NULL
best_sig <- -Inf
best_aic <- Inf
if(qq == 0){
n <- 0
} else{
n <- 1
}
for (i in 1:pp) {
for (j in n:qq) {
# GARCH model
formula <- as.formula(paste("~ garch(", i, ",", j, ")", sep = ""))
archSpec <- create_spec(i, j, arma.order, mean=mean)
archFit <- ugarchfit(spec = archSpec, data = data, solver = "hybrid")
if (anyNA(coef(archFit)) || anyNA(archFit@fit[["tval"]])) {
break
}
aic <- infocriteria(archFit)[[1]] # AIC
p_values <- archFit@fit$robust.matcoef[, 4] # P-value
# persentase banyaknya param yang signif
sig <- sum(p_values < 0.05) / length(p_values)
sig_param <- names(p_values[p_values < 0.05])
non_sig_param <- names(p_values[p_values >= 0.05])
model <- if (qq == 0) {
paste0("ARCH(", i, ")")
} else {
paste0("GARCH(", i, ",", j, ")")
}
arch.model <- rbind(
arch.model,
data.frame(
Model = model,
Signif = sig,
AIC = aic,
Sig.par = paste(sig_param, collapse = ", "),
Non_Sig.par = paste(non_sig_param, collapse = ", ")
)
)
# Identifikasi best model
if (sig > best_sig) {
if (aic < best_aic) {
best_model <- model
best_sig <- sig
best_aic <- aic
}
}
}
}
return(list(arch.model = arch.model, best_model = best_model))
}
Model ARCH tanpa Intersep
model.arch1 <- arch.garch(data = train_nvda.diff, pp = 12, qq = 0, arma.order = arma.order)
datatable(model.arch1$arch.model, filter = 'top',
options = list(pageLength = 6))
cat("Best Model:", model.arch1$best_model)
## Best Model: ARCH(1)
Setelah mencoba iterasi lag 1 sampai lag 12, terlihat bahwa hanya ada satu model yang semua parameternya signifikan yakni ARCH(1). Karena itu tidak perlu membandingkan nilai AIC. Sehingga model terbaik adalah ARIMA(0,1,1)+ARCH(1).
Model ARCH dengan intersep
model.arch2 <- arch.garch(data=train_nvda.diff, pp=12, qq=0, arma.order=arma.order,
mean=TRUE)
datatable(model.arch2$arch.model, filter = 'top',
options = list(pageLength = 6))
cat("Best Model:", model.arch2$best_model)
## Best Model: ARCH(1)
Sama seperti model ARCH tanpa intersep, model terbaiknya adalah ARCH(1).
#Model dengan intersep
archSpec <- ugarchspec(
variance.model=list(model="sGARCH", garchOrder=c(1,0)),
mean.model=list(armaOrder = arma.order, include.mean = TRUE) )
model.arch_nvda1 <- ugarchfit(spec=archSpec, data=train_nvda.diff, solver="hybrid")
coef(model.arch_nvda1)
## mu ma1 omega alpha1
## -5.087292e-05 4.308273e-02 7.454659e-07 2.989799e-01
#Model tanpa intersep
archSpec <- ugarchspec(
variance.model=list(model="sGARCH", garchOrder=c(1,0)),
mean.model=list(armaOrder = arma.order, include.mean = FALSE) )
model.arch_nvda2 <- ugarchfit(spec=archSpec, data=train_nvda.diff, solver="hybrid")
coef(model.arch_nvda2)
## ma1 omega alpha1
## 4.514680e-02 7.495716e-07 2.961780e-01
Model yang dihasilkan adalah seperti diatas.
6.1.3 Plot
plot.ts(ts(model.arch_nvda1@fit$residuals), col=col.nvda[1], lty=1,
xlab="Periode (Hari)", ylab="Conditional Variance",
lwd = 2, main="ARIMA(0,1,1)-ARCH(1) NVDA", cex.main = 1.5)
Terlihat bahwa ada periode-periode yang volatilitas nya tinggi.
6.2 GARCH
Referensi : GARCH models with R programming : a practical example with TESLA stock
Ya walaupun diatas sudah dijelaskan “tidak usah melanjutkan ke garch”, tapi kami hanya ingin mencoba. “Apa yang terjadi jika kita tetap melakukan model GARCH?”. Model terbaik dari GARCH(p,q) ditentukan dari banyaknya paramter yang signifikan dan AIC nya yang terkecil.
model.garch <- arch.garch(data = train_nvda.diff, pp = 5, qq = 5, arma.order = arma.order)
datatable(model.garch$arch.model, filter = 'top',
options = list(pageLength = 6))
cat("Best Model:", model.garch$best_model)
## Best Model: GARCH(1,1)
Terlihat bahwa memang data tidak cocok untuk model GARCH(p,q) karena dari rentang lag 1 sampai 12 yang dihasilkan model terbaiknya adalah GARCH(12,12). Namun jika pada modelnya, hanya persentase signifikannya sangat kecil. Sehingga tidak mungin untuk dijadikan model terbaik.
6.3 Model Terbaik
Sehingga diperoleh model terbaik yakni ARIMA(0,1,1)+ARCH(1).
plot(model.arch_nvda1, which="all")
##
## please wait...calculating quantiles...
Diperoleh mean model yaitu ARIMA(0,1,1) dan varian model yaitu ARCH(1) untuk digunakan dalam peralaman, pada uji dignostik model pada plot QQ-norm terlihat bahwa sebaran residual tidak mengikuti garis lurus sehingga dapat disimpulkan bahwa residual tidak menyebar normal.
Box.test((sisaan_nvda.da)^2, type = "Ljung")
##
## Box-Ljung test
##
## data: (sisaan_nvda.da)^2
## X-squared = 4.2891, df = 1, p-value = 0.03836
7 Akurasi
Akan dilakukan peramalan sebanyak data uji untuk mengukur akurasi.
Peramalan dilakukan menggunakan fungsi forecast() . Contoh
peramalan berikut ini dilakukan untuk \(h\) hari ke depan. Dengan \(h\) adalah banyaknya data test per
perusahaan.
#--Forecast--
#ARIMA
ramalan_nvda.da <- forecast::forecast(model_nvda.da, h = nrow(test_nvda))
data.ramalan_nvda.da <- ramalan_nvda.da$mean
#ARCH + Intersep
ramalan_nvda.arch1 = ugarchforecast(model.arch_nvda1, data = train_nvda.diff,
n.ahead = nrow(test_nvda))
data.ramal_nvda.arch1 <- ramalan_nvda.arch1@forecast[["seriesFor"]]
#ARCH - Intersep
ramalan_nvda.arch2 = ugarchforecast(model.arch_nvda2, data = train_nvda.diff,
n.ahead = nrow(test_nvda))
data.ramal_nvda.arch2 <- ramalan_nvda.arch2@forecast[["seriesFor"]]
7.1 Plot Ramalan
Model ARIMA(0,1,1)
#Plot
plot(ramalan_nvda.da, main="NVDA", xlab="Periode(Hari)", cex.main=4,
lwd=2.5, col = col.nvda[1])
Berdasarkan hasil plot ramalan di atas, dapat dilihat bahwa semua ramalan cenderung stabil hingga akhir periode.
Model ARIMA(0,1,1)+ARCH(1)
plot(ramalan_nvda.arch1, which= 1)
Berdasarkan hasil plot ramalan di atas, dapat dilihat bahwa semua ramalan cenderung stabil hingga akhir periode.
#--Diff inv--
#ARIMA
ramal_nvda <- diffinv(data.ramalan_nvda.da, differences=1) +
train_nvda.ts.new[length(train_nvda.ts)] #nilai akhir data latih
#ARCH + Intersep
ramal_nvda.arch1 <- diffinv(data.ramal_nvda.arch1, differences=1) +
train_nvda.ts.new[length(train_nvda.ts)] #nilai akhir data latih
#ARCH - Intersep
ramal_nvda.arch2 <- diffinv(data.ramal_nvda.arch2, differences=1) +
train_nvda.ts.new[length(train_nvda.ts)] #nilai akhir data latih
#--Transformasi Balik--
#ARIMA
ramal_nvda <- 1/(ramal_nvda)^2
#ARCH + Intersep
ramal_nvda.arch1 <- 1/(ramal_nvda.arch1)^2
#ARCH - Intersep
ramal_nvda.arch2 <- 1/(ramal_nvda.arch2)^2
#Data Frame
ramal <- data.frame(
Date = data2 %>% dplyr::select(Date) ,
`ARIMA` = c(as.vector(train_nvda.ts), ramal_nvda[-1] ),
`ARCH+Intersep` = c(as.vector(train_nvda.ts), ramal_nvda.arch1[-1] ),
`ARCH-Intersep` = c(as.vector(train_nvda.ts), ramal_nvda.arch2[-1] ),
Name = rep(c("NVDA"),
each = nrow(train_nvda) + nrow(test_nvda) )
)
colnames(ramal) <- c("Date", "ARIMA", "ARCH+Intersep", "ARCH-Intersep", "Name")
chart <-
ggplot() +
# Akhir data
geom_vline( #Buat garis batas data
xintercept = as.numeric(as.Date( max(data2$Date) )),
linetype = "dotted", color = "red") +
#Label Data Asli
annotate( "rect", alpha=.1, fill="#4EC2C1",
xmin=as.Date(min(data2$Date)),
xmax=as.Date(data2$Date[nrow(train_nvda)]),
ymin=max(nvda$`Adj Close`) * 1.25, ymax=Inf ) +
annotate( "text", color="#4EC2C1",
x = as.Date(data2$Date[ceiling(nrow(train_nvda)/2)]),
y = max(nvda$`Adj Close`) * 1.3,
label = "Data Asli", size=10) +
annotate( "text", color="#4EC2C1",
x = as.Date(data2$Date[ceiling(nrow(train_nvda)/2)]),
y = max(nvda$`Adj Close`) * 1.2,
label = paste0(nrow(train_nvda)," Hari (",
round(nrow(train_nvda)/nrow(nvda)*100, 1),
"%)"), size=10) +
#Label Data Ramal
annotate( "rect", alpha=.1, fill="violetred",
xmin=as.Date(data2$Date[nrow(train_nvda)]),
xmax=as.Date(max(data2$Date)),
ymin=max(nvda$`Adj Close`) * 1.25, ymax=Inf ) +
annotate( "text", color="violetred",
x = as.Date(data2$Date[ceiling( .89*length(data2$Date)/6)]) ,
y = max(nvda$`Adj Close`) * 1.3,
label = "Data Ramal", size=10) +
annotate( "text", color="violetred",
x = as.Date(data2$Date[ceiling( .89*length(data2$Date)/6)]) ,
y = max(nvda$`Adj Close`) * 1.2,
label = paste0(nrow(test_nvda)," Hari (",
round(nrow(test_nvda)/nrow(nvda)*100, 1),
"%)"), size=10) +
geom_vline( #Buat garis batas data
xintercept = as.Date(data2$Date[nrow(train_nvda)]) ,
linetype = "dotted", color = "black", linewidth = 2) +
#Time Series
#Data RAMAL ARIMA
#Dari data asli & ramal
geom_line(data = ramal, linewidth=2,
aes(x = Date, y = `ARIMA`, col = "ARIMA")) +
#Data Ramal ARCH + Intersep
#Dari data ramal saja
geom_line(data = ramal[nrow(train_nvda):nrow(nvda),], linewidth=2,
aes(x = Date, y = `ARCH+Intersep`, col = "ARCH+")) +
#Data Ramal ARCH - Intersep
#Dari data ramal saja
geom_line(data = ramal[nrow(train_nvda):nrow(nvda),], linewidth=2,
aes(x = Date, y = `ARCH-Intersep`, col = "ARCH-")) +
#Data Test
geom_line(data = test_nvda, linewidth=2,
aes(x = Date, y = `Adj Close`, col = "Test")) +
#Penanda
geom_point(data = tail(train_nvda, 1), alpha = .5,
aes(x = Date, y = `Adj Close`), stroke=2,
size = 15, shape = 21, color = "black", fill="violetred") +
scale_colour_manual(values = c("ARIMA"=col.nvda[1],"ARCH+"="#539282",
"ARCH-"="#93C6B9","Test"=col.nvda[2]
)) +
theme.ts + #THeme
labs(x = "\nPeriode (Tahun)", y='Harga Saham (USD)',
title = "Akurasi Ramalan Saham NVDA",
subtitle = "Perbandingan Hasil Ramalan dengan data Asli\n") +
#Axis
coord_cartesian(clip = "off"
) +
scale_x_date( #Sumbu x
date_breaks = "1 year", # Menampilkan label setiap tahun
date_labels = "%Y", # Format label tahun
limits = c(as.Date(min(ramal$Date)),
as.Date(max(ramal$Date)) + 80)
) +
scale_y_continuous( #Sumbu y
labels = scales::dollar_format(prefix = "$") #tambahin dolar
) +
#Bagi Data Train & Test
geom_segment(aes(x = as.Date(max(train_nvda$Date)),
xend = as.Date(max(train_nvda$Date)),
y = max(nvda$`Adj Close`)* .58,
yend = max(nvda$`Adj Close`)*.7),
arrow = arrow(type = "closed", length = unit(0.1, "inches")),
lineend = "round", color = "gray60", size=1.5) +
geom_richtext(
data = data.frame(x = as.Date(max(train_nvda$Date)),
y = max(nvda$`Adj Close`)*.7,
label = max(train_nvda$Date) ),
aes(x, y, label = label), size = 8, color = "white",
fill = "gray60", box.color = "white", parse = TRUE
) +
#Data Uji
geom_segment(aes(x = as.Date(max(data2$Date)),
xend = as.Date(max(data2$Date))+60,
y = nvda$`Adj Close`[nrow(nvda)],
yend = nvda$`Adj Close`[nrow(nvda)]),
arrow = arrow(type = "closed", length = unit(.1, "inches")),
lineend = "round", color = "gray60", size=1.5) +
geom_richtext(
data = data.frame(x = as.Date(max(data2$Date))+60,
y = nvda$`Adj Close`[nrow(nvda)],
label = "Data Asli" ),
aes(x, y, label = label), size = 8, color = "white",
fill = "gray60", box.color = "white", parse = TRUE
) +
# Data Ramal ARIMA
geom_segment(aes(x = as.Date(max(data2$Date)),
xend = as.Date(max(data2$Date)) + 60,
y = max(ramal$`ARIMA`),
yend = max(ramal$`ARIMA`)),
arrow = arrow(type = "closed", length = unit(.1, "inches")),
lineend = "round", color = "#D02A49", size=1.5) +
geom_richtext(
data = data.frame(x = as.Date(max(data2$Date)) +60,
y = max(ramal$`ARIMA`),
label = "ARIMA" ),
aes(x, y, label = label), size = 8, color = "white",
fill = "#D02A49", box.color = "white", parse = TRUE
) +
# Data Ramal ARCH + Intersep
geom_segment(aes(x = as.Date(max(data2$Date)),
xend = as.Date(max(data2$Date)) + 70,
y = max(ramal$`ARCH+Intersep`),
yend = max(ramal$`ARCH+Intersep`)),
arrow = arrow(type = "closed", length = unit(.1, "inches")),
lineend = "round", color = "#4B75BA", size=1.5) +
geom_richtext(
data = data.frame(x = as.Date(max(data2$Date)) +70,
y = max(ramal$`ARCH+Intersep`),
label = "ARCH dengan <br> Intersep" ),
aes(x, y, label = label), size = 8, color = "white",
fill = "#4B75BA", box.color = "white", parse = TRUE
) +
# Data Ramal ARCH - Intersep
geom_segment(aes(x = as.Date(max(data2$Date)),
xend = as.Date(max(data2$Date)) + 70,
y = ramal$`ARCH-Intersep`[nrow(nvda)],
yend = ramal$`ARCH-Intersep`[nrow(nvda)]),
arrow = arrow(type = "closed", length = unit(.1, "inches")),
lineend = "round", color = "#DC4F26", size=1.5) +
geom_richtext(
data = data.frame(x = as.Date(max(data2$Date)) +70,
y = ramal$`ARCH-Intersep`[nrow(nvda)],
label = "ARCH tanpa <br> Intersep" ),
aes(x, y, label = label), size = 8, color = "white",
fill = "#DC4F26", box.color = "white", parse = TRUE
) +
#Tanggal akhir
geom_segment(aes(x = as.Date(max(data2$Date)),
xend = as.Date(max(data2$Date)),
y = -Inf,
yend = max(nvda$`Adj Close`)*.15),
arrow = arrow(type = "closed", length = unit(.1, "inches")),
lineend = "round", color = "gray60", size=1.5) +
geom_richtext(
data = data.frame(x = as.Date(max(data2$Date)),
y = max(nvda$`Adj Close`)*.15,
label = max(data2$Date) ),
aes(x, y, label = label), size = 8, color = "white",
fill = "gray60", box.color = "white", parse = TRUE
)
chart
#Export Chart
ggsave("06_Akurasi_model.png", chart, path = export.chart,
dpi = 300, height = 12, width = 27)
Dapat dilihat bahwa rata-rata harga saham Amazon diramalkan akan cenderung sedikit menurun setiap periodenya. Di sisi lain, saham NVIDIA diramalkan akan menaik pesat setiap periodenya. Sedangkan harga saham lainnya cenderung naik sedikit.
7.2 Perbandingan Data Aktual dengan Data Ramal
#Perbandingan
perbandingan_nvda.da <- matrix(data=c(head(test_nvda.ts, n=nrow(test_nvda)),
ramal_nvda[-1],
ramal_nvda.arch1[-1],
ramal_nvda.arch2[-1]),
nrow = nrow(test_nvda), ncol = 4)
colnames(perbandingan_nvda.da) <- c("Aktual","ARIMA",
"ARCH+Intersep","ARCH-Intersep")
datatable(perbandingan_nvda.da, filter = 'top',
options = list(pageLength = 5))
7.3 Akurasi
#Akurasi Model ARIMA
akurasi.nvda <- accuracy(ts(ramal_nvda[-1]),
head(test_nvda.ts, n=nrow(test_nvda) )) %>%
as.data.frame() %>%
cbind(Model = ramalan_nvda.da[["method"]]) %>%
relocate(Model, .before = 1)
row.names(akurasi.nvda) <- "ARIMA"
#Akurasi Model ARCH + Intersep
akurasi.nvda.arch1 <- accuracy(ts(ramal_nvda.arch1[-1]),
head(test_nvda.ts, n=nrow(test_nvda) )) %>%
as.data.frame() %>%
cbind(Model = model.arch1$best_model) %>%
relocate(Model, .before = 1)
row.names(akurasi.nvda.arch1) <- "ARCH+Intersep"
#Akurasi Model ARCH - Intersep
akurasi.nvda.arch2 <- accuracy(ts(ramal_nvda.arch2[-1]),
head(test_nvda.ts, n=nrow(test_nvda) )) %>%
as.data.frame() %>%
cbind(Model = model.arch2$best_model) %>%
relocate(Model, .before = 1)
row.names(akurasi.nvda.arch2) <- "ARCH-Intersep"
akurasi <- rbind(akurasi.nvda, akurasi.nvda.arch1, akurasi.nvda.arch2) %>%
mutate(Keterangan = case_when(
MAPE < 10 ~ "Sangat Baik",
MAPE >= 10 & MAPE <= 20 ~ "Baik",
MAPE > 20 & MAPE <= 50 ~ "Layak",
MAPE > 50 ~ "Tidak Akurat"
)) %>% relocate(Keterangan, Model, MAPE)
datatable(akurasi, filter = 'top',
options = list(pageLength = 7))
Dari tabel diatas terlihat bahwa model ARIMA(0,1,1)-ARCH(1) dengan Intersep berkatergori “Baik” dengan nilai MAPE nya sebesar \(15.64\%\) (kurang dari \(20\%\)).
8 Peramalan
8.1 Stasioneritas
8.1.1 Ragam
index <- seq(1:nrow(nvda)) #Beda titik potong
bc =boxcox(nvda.ts~index, lambda=seq(-2, 4, by=.01))
title("NVDA", cex.main = 2)
lambda_nvda.ramal <- bc$x[which.max(bc$y)] #Nilai Rounded Lambda
sk_nvda.ramal <- bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)] #SK
lsk.ramal <- data.frame(
Lambda = c(lambda_nvda.ramal ),
Batas.Bawah = c( min(sk_nvda.ramal) ),
Batas.Atas = c(max(sk_nvda.ramal) )
) %>%
mutate(across(where(is.numeric), ~round(., 2)))
lsk.ramal$Keterangan <- apply(lsk.ramal, 1, function(row) {
ifelse(row["Batas.Bawah"] <= 1 && row["Batas.Atas"] >= 1,
"Ragam Stasioner", "Ragam Tidak Stasioner")
})
rownames(lsk.ramal) <- c("NVDA")
datatable(lsk.ramal, filter = "none")
Data tidak Stasioner dalam ragam.
tran.ramal <- lsk.ramal %>%
mutate(
Transformasi = trans(Lambda),
`Transformasi Balik` = trans_balik(Lambda)
) %>% dplyr::select(Lambda, Transformasi, `Transformasi Balik`)
datatable(tran.ramal, filter = 'none')
Berikut transformasi yang cocok.
8.1.2 Transformasi
#Transformasi
nvda.ts.new <- 1/sqrt(nvda.ts)
#Plot BoxCox
index <- seq(1:nrow(nvda))
bc <-boxcox(nvda.ts.new ~index, lambda=seq(-2, 4, by=.01))
title("NVDA", cex.main=2)
lambda_nvda.ramal.new <- bc$x[which.max(bc$y)] #Nilai Rounded Lambda
sk_nvda.ramal.new <- bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)] #SK
lsk.ramal.new <- data.frame(
Lambda = c(lambda_nvda.ramal.new ),
Batas.Bawah = c( min(sk_nvda.ramal.new) ),
Batas.Atas = c( max(sk_nvda.ramal.new) )
) %>%
mutate(across(where(is.numeric), ~round(., 2)))
lsk.ramal.new$Keterangan <- apply(lsk.ramal.new, 1, function(row) {
ifelse(row["Batas.Bawah"] <= 1 && row["Batas.Atas"] >= 1,
"Ragam Stasioner", "Ragam Tidak Stasioner")
})
rownames(lsk.ramal.new) <- c("NVDA")
datatable(lsk.ramal.new, filter = 'none')
Terlihat bahwa data sudah stasioner dalam ragam.
8.1.3 Rataan & Differencing
#Diff
nvda.diff <- diff(nvda.ts.new, differences = 1)
tseries::adf.test(nvda.diff)
##
## Augmented Dickey-Fuller Test
##
## data: nvda.diff
## Dickey-Fuller = -8.9643, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary
if (tseries::adf.test(nvda.diff)[["p.value"]] < 0.05) {
cat("Karena p-value < 0.05, Maka Rataan Stasioner")
} else {
cat("Karena p-value > 0.05, Maka Rataan Tidak Stasioner")
}
## Karena p-value < 0.05, Maka Rataan Stasioner
Setelah di differencing, data sudah stasioner.
8.2 Model
#Model
archSpec <- ugarchspec(
variance.model=list(model="sGARCH", garchOrder=c(1,0)), #ARCH(1)
mean.model=list(armaOrder = arma.order, include.mean = TRUE) ) #ARIMA(0,1,1) + Intersep
model.ramal <- ugarchfit(spec=archSpec, data=nvda.diff, solver="hybrid")
#--Forecast--
nDate <- seq(from = max(data2$Date) + 1,
to = as.Date("2023-12-31"), #Dampe akhir tahun
by = "days")
ramalan_nvda <- ugarchforecast(model.ramal, data = nvda.diff,
n.ahead = length(nDate))
data.ramal_nvda <- ramalan_nvda@forecast[["seriesFor"]]
data.ramal_nvda.bb <- ramalan_nvda@forecast[["seriesFor"]]/
-( ramalan_nvda@forecast[["sigmaFor"]] )^2
data.ramal_nvda.ba <- ramalan_nvda@forecast[["seriesFor"]]/
( ramalan_nvda@forecast[["sigmaFor"]] )^2
#--Diff inv--
#Mean
ramal_nvda <- diffinv(data.ramal_nvda, differences=1) +
nvda.ts.new[length(nvda.ts)] #nilai akhir data latih
#BB
ramal_nvda.bb <- diffinv(data.ramal_nvda.bb, differences=1) +
nvda.ts.new[length(nvda.ts)] #nilai akhir data latih
#BA
ramal_nvda.ba <- diffinv(data.ramal_nvda.ba, differences=1) +
nvda.ts.new[length(nvda.ts)] #nilai akhir data latih
#--Transformasi Balik--
ramal_nvda <- 1/(ramal_nvda)^2
ramal_nvda.bb <- 1/(ramal_nvda.bb)^2
ramal_nvda.ba <- 1/(ramal_nvda.ba)^2
#Data Frame
data_ramal <- data.frame(
Date = seq(from = min(data2$Date),
to = max(nDate), by="days") ,
`Adj Close` = c(as.vector(nvda.ts), ramal_nvda[-1] ),
`BB` = c(as.vector(nvda.ts), ramal_nvda.bb[-1] ),
`BA` = c(as.vector(nvda.ts), ramal_nvda.ba[-1] ),
Name = rep(c("NVDA"),
each = length(nvda.ts) + length(nDate) )
)
colnames(data_ramal) <- c("Date", "Adj Close", "BB", "BA", "Name")
8.3 Plot
harga_akhir <- nvda$`Adj Close`[nrow(nvda)]
harga_ramal <- data_ramal$`Adj Close`[nrow(data_ramal)]
perc <- (harga_ramal-harga_akhir)/harga_akhir * 100
chart <-
ggplot() +
# Akhir data
geom_vline( #Buat garis batas data
xintercept = as.numeric(as.Date( max(data_ramal$Date) )),
linetype = "dotted", color = "red") +
#Label Data Asli
annotate( "rect", alpha=.1, fill="#4EC2C1",
xmin=as.Date("2022-06-01"),
xmax=as.Date(nvda$Date[nrow(nvda)]),
ymin=max(nvda$`Adj Close`) * 1.25, ymax=Inf ) +
annotate( "text", color="#4EC2C1",
x = as.Date(data_ramal$Date[ceiling(nrow(data_ramal)*.7)]),
y = max(nvda$`Adj Close`) * 1.3,
label = "Saham Historis", size=10) +
annotate( "text", color="#4EC2C1",
x = as.Date(data_ramal$Date[ceiling(nrow(data_ramal)*.7)]),
y = max(nvda$`Adj Close`) * 1.2,
label = paste0(" "), size=10) +
#Label Data Ramal
annotate( "rect", alpha=.1, fill="violetred",
xmin=as.Date(data_ramal$Date[nrow(nvda)]),
xmax=as.Date(max(data_ramal$Date)),
ymin=-Inf, ymax=Inf ) +
annotate( "text", color="violetred",
x = as.Date(max(nDate)-25) ,
y = min(data2$`Adj Close`) * 10,
label = "Data\nRamal", size=10) +
geom_vline( #Buat garis batas data
xintercept = as.Date(data_ramal$Date[nrow(nvda)]) ,
linetype = "dotted", color = "black", linewidth = 2) +
#Time Series
#Data RAMAL
geom_line(data = data_ramal, linewidth=2,
aes(x = Date, y = `Adj Close`, col = "NVDA")) +
#Penanda
geom_point(data = tail(nvda, 1), alpha = .5,
aes(x = Date, y = `Adj Close`), stroke=2,
size = 15, shape = 21, color = "black", fill="violetred") +
scale_colour_manual(values = cols) +
theme.ts + #THeme
labs(x = "\nPeriode (Tahun)", y='Harga Saham (USD)',
title = "Ramalan Saham NVDA",
subtitle = "Peramalan Saham NVDA Hingga Akhir Tahun 2023\n") +
#Axis
coord_cartesian(clip = "off"
) +
scale_x_date( #Sumbu x
date_breaks = "1 year", # Menampilkan label setiap tahun
date_labels = "%Y", # Format label tahun
limits = c(as.Date("2022-06-01"),
as.Date(max(data_ramal$Date)) + 40)
) +
scale_y_continuous( #Sumbu y
labels = scales::dollar_format(prefix = "$") #tambahin dolar
) +
#Tanggal akhir data asli
geom_richtext(
data = data.frame(x = as.Date(max(nvda$Date)),
y = max(nvda$`Adj Close`)*.75,
label = max(nvda$Date) ),
aes(x, y, label = label), size = 8, color = "white",
fill = "gray60", box.color = "white", parse = TRUE
) +
# Harga akhir Data Asli
geom_segment(aes(x = as.Date(max(nvda$Date)),
xend = as.Date(max(nvda$Date))- 30,
y = nvda$`Adj Close`[nrow(nvda)],
yend = max(data_ramal$`Adj Close`)*.9 ) ,
arrow = arrow(type = "closed", length = unit(.1, "inches")),
lineend = "round", color = "#4B75BA", size=1.5) +
geom_richtext(
data = data.frame(x = as.Date(max(nvda$Date)) - 30,
y = max(data_ramal$`Adj Close`) *.9,
label = paste0("$ ", harga_akhir %>% round(.,2)) ),
aes(x, y, label = label), size = 8, color = "white",
fill = "#4B75BA", box.color = "white", parse = TRUE
) +
# Harga akhir Data Ramal
geom_segment(aes(x = as.Date(max(data_ramal$Date)),
xend = as.Date(max(data_ramal$Date)) + 30,
y = max(data_ramal$`Adj Close`),
yend = max(data_ramal$`Adj Close`)),
arrow = arrow(type = "closed", length = unit(.1, "inches")),
lineend = "round", color = "#4B75BA", size=1.5) +
geom_richtext(
data = data.frame(x = as.Date(max(data_ramal$Date)) +30,
y = max(data_ramal$`Adj Close`),
label = paste0("$ ", harga_ramal %>% round(.,2)) ),
aes(x, y, label = label), size = 8, color = "white",
fill = "#4B75BA", box.color = "white", parse = TRUE
) +
# Harga akhir data ramal beda brp persen
annotate( "text", color=ifelse(perc>0, "green4", "red3"),
x = as.Date(max(data_ramal$Date)) +30,
y = max(data_ramal$`Adj Close`) *.94,
label = paste0(ifelse(perc>0, "+ ", "- "), perc %>% round(.,2), "%"),
size=8 ) +
#Tanggal akhir ramal
geom_segment(aes(x = as.Date(max(data_ramal$Date)),
xend = as.Date(max(data_ramal$Date)),
y = -Inf,
yend = max(nvda$`Adj Close`)*.15),
arrow = arrow(type = "closed", length = unit(.1, "inches")),
lineend = "round", color = "gray60", size=1.5) +
geom_richtext(
data = data.frame(x = as.Date(max(data_ramal$Date)),
y = max(nvda$`Adj Close`)*.15,
label = max(data_ramal$Date) ),
aes(x, y, label = label), size = 8, color = "white",
fill = "gray60", box.color = "white", parse = TRUE
)
chart
#Export Chart
ggsave("07_Ramalan.png", chart, path = export.chart,
dpi = 300, height = 12, width = 27)
Model ARIMA(0,1,1)-ARCH(1) dengan intersep digunakan untuk meramalkan hingga periode akhir tahun 2023. Pada Gambar 12 terlihat bahwa adjusted close saham NVIDIA diramalkan akan menyentuh angka seharga 555,95 USD atau naik sebesar \(15,02\%\) hingga akhir tahun 2023.
9 Penutup
#Export Data
install_load('openxlsx')
#Model Tentatif
write.xlsx(list("Tentatif" = model.tentaif_nvda,
"Drift1" = drift1_nvda,
"Drift2" = drift2_nvda,
"Drift3" = drift3_nvda,
"Drift4" = drift4_nvda,
"Arch-Intersep" = model.arch1$arch.model,
"Arch+Intersep" = model.arch2$arch.model,
"Garch" = model.garch$arch.model
),
file = "Model_NVDA.xlsx")
9.1 Kesimpulan
Kesimpulan yang dapat ditarik pada penelitian ini adalah model ARIMA(0,1,1)-ARCH(1) dengan intersep merupakan model terbaik yang dapat meramalkan adjusted close saham NVIDIA dalam kurung waktu 11 November 2023 hingga 31 Desember 2023. Nilai MAPE yang didapatkan dari model terbaik adalah sebesar \(15.64\%\). Adjusted close saham NVIDIA diramalkan akan naik hingga menyentuh angka \(555.96\) USD atau naik sebesar \(15.02\%\). Hal ini menunjukkan bahwa jika ingin membeli saham harian NVIDIA pada kurung waktu 11 November 2023 hingga 31 Desember 2023 merupakan hal yang tepat dikarenakan saham akan naik.
9.2 Referensi
Forecasting: Principles and Practice, Rob J Hyndman and George Athanasopoulos. Monash University, Australia.
BOOSTEDML : Time series, Articles on Statistics and Machine Learning for Healthcare.